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Abstract 

A  one-dimensional  rigorous  process  model  of  a  single-cell  direct  methanol  fuel  cell  (DMFC)  is  presented.  Multi-component  mass  transport 
in  the  diffusion  layers  and  the  polymer  electrolyte  membrane  (PEM)  is  described  using  the  generalised  Maxwell-Stefan  (MS)  equation  for 
porous  structures.  In  the  PEM,  also  local  swelling  behaviour  and  non-idealities  are  accounted  for  by  a  Flory-Huggins  model  for  the  activities 
of  the  mobile  species  inside  the  pores  of  the  PEM.  Phase  equilibria  between  the  pore  liquid  inside  the  PEM  and  those  inside  the  pores  of  both 
catalyst  layer  are  formulated  based  on  literature  data  and  activity  models.  Although  two-phase  behaviour  in  both  diffusion  layers  is  neglected, 
the  model  shows  good  agreement  to  own  experimental  data  over  a  wide  range  of  operating  conditions,  with  respect  to  methanol  and  water 
crossover  fluxes  as  well  as  to  current-voltage  characteristics.  Only  for  very  low  current  densities  and  in  the  limiting  current  regime  significant 
deviations  between  model  and  experiments  are  found. 
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Abbreviations:  A,  anode  compartment  (supply  channel  structure);  AC, 
anode  catalyst  layer;  ACP,  polymer-phase  within  (AC);  AD,  anode  diffu¬ 
sion  layer;  BV,  Butler- Volmer  type  rate  expression  (Table  1);  C,  cathode 
compartment  (supply  channel  structure);  CC,  cathode  catalyst  layer;  CCP, 
polymer-phase  within  (CC);  CD,  cathode  diffusion  layer;  dc,  drag  coeffi¬ 
cient  mass  transport  model  (Table  1);  dyn.,  dynamic  (Table  1);  DMFC,  direct 
methanol  fuel  cell;  eff.,  effective  (Table  1);  F,  Fick  diffusion  model  (Table  1); 
g,  gas;  irrev.,  irreversible  (Table  1);  1,  liquid;  M,  membrane  (PEM);  MEA, 
membrane  electrode  assembly  (DMFC  core  component);  MS,  Maxwell- 
Stefan  mass  transport  model  (Table  1);  NP,  Nernst-Planck  mass  transport 
model  (Table  1);  PEM,  polymer  electrolyte  membrane;  PEMFC,  polymer 
electrolyte  membrane  fuel  cell;  PTFE,  polytetrafluoeethylene,  TEFLON™; 
scbm,  standard  cubic  metre  (m3  ideal  gas  at  T  =  25°C,  p  =  1  bar);  S, 
Schlogl  approach  to  convective  mass  transport  (Table  1);  SD,  surface  diffu¬ 
sion  mass  transport  model  (Table  1);  s.s.,  steady-state  (Table  1);  TD,  ther¬ 
modynamics  (Table  1 ) 
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1.  Introduction 

Since  first  direct  methanol  fuel  cell  (DMFC)  systems  are 
commercially  available  for  special  outdoor  and  professional 
applications  (e.g.  from  SmartFuelCell  GmbH,  Germany),  the 
prospects  of  this  type  of  fuel  cell  are  gaining  an  even  higher 
interest  in  the  consumer  electronics  industry.  The  first  com¬ 
mercial  systems  prove  the  practical  value  and  applicability 
of  the  DMFC,  but  also  the  challenges  still  ahead.  The  fuel 
efficiencies  and  the  power  densities  are  still  very  low,  and 
a  reliable  dynamic  operation  without  a  significant  buffer  for 
electrical  energy  (like  e.g.  a  battery  or  a  supercapacitor)  has 
not  yet  been  reported. 

The  necessary  improvements  in  DMFC  performance  and 
operation  do  not  only  place  the  demand  for  better  materi¬ 
als,  i.e.  better  catalysts  and  a  polymer  electrolyte  membrane 
(PEM)  less  or  even  impermeable  for  methanol  and  water, 
but  also  for  sophisticated  controller  designs.  Only  the  lat- 


0378-7753/$  -  see  front  matter  ©  2005  Elsevier  B.V.  All  rights  reserved, 
doi:  10. 1016/j.jpowsour.2005.02.036 


436 


T.  Schultz,  K.  Sundmacher  /  Journal  of  Power  Sources  145  (2005)  435-462 


Nomenclature 

a 

activity  (-) 

au2o 

water  vapour  activity  (-) 

As~ 

cell  cross-sectional  area  (m2) 

At 

parameter  in  empirical  correlations  (-) 

B 

transport  matrix 

Bo 

permeability  coefficient  (mm2) 

Bi 

parameter  in  empirical  correlations  (-) 

c 

molar  concentration  in  fluid-phase  (molm-3) 

c 

molar  pseudo-concentration  w.r.t.  total  volume 
(in porous  structures  only)  (molm-3) 

Ci 

parameter  in  empirical  correlations  (-) 

cP 

mass-based  heat  capacity  at  constant  pressure 
(Jkg-1K-1) 

Cp 

molar  heat  capacity  at  constant  pressure 
(Jmol-1K-1) 

ArCp 

molar  heat  capacity  change  of  reaction  at  con¬ 
stant  pressure  (J  mol-1K-1) 

d 

thickness,  diameter  (m) 

D 

diffusion  coefficient  (mm2  s-1) 

D 

Maxwell-Stefan  binary  diffusion  coefficient 
(m2  s-1) 

D, 

parameter  in  empirical  correlations  (-) 

e 

enthalpy  flux  density  (J  m-2  s-1) 

Ea 

activation  energy  (J  mol-1) 

Ei 

parameter  in  empirical  correlations  (-) 

F 

Faraday’s  constant,  F  — 

96485  As  mol- '(As  mol-1) 

ApG 

Gibbs  energy  of  formation  (from  the  elements) 

(J  mol-1) 

ArG 

Gibbs  energy  of  reaction  (J  mol-1) 

h 

specific  enthalpy  (J  mol-1) 

AF  H 

enthalpy  of  formation  (from  the  elements) 

(J  mol-1) 

Ar  H 

reaction  enthalpy  (J  mol-1) 

i 

current  density  (A  m-2) 

j 

individual  molar  flux  density  (mol  m-2  s-1) 

k 

index  for  control  volumes  (discretised  model) 

(-) 

Li 

friction  terms  (s  m-2) 

m 

mass  flux  density  (kg  m-2  s-1) 

M 

mass  (kg) 

M 

molar  mass  (kg  mol-1) 

n 

overall  molar  flux  density  (mol  m-2  s-1) 

N 

number  of  moles  (mol) 

N 

mole  density  (loading,  used  only  in  polymer 
material)  (mol  m-2) 

Nm.cu 

number  of  chain  units  between  two  polymer 
cross-links  (-) 

p 

pressure  (Pa) 

Psat 

saturation  pressure  (Pa) 

P 

parachor  (cm3  g0  25  s-0  5) 

q 

heat  flux  density  (due  to  thermal  conduction) 

(J  m-2  s-1) 

Q 

charge  (C  =  As) 

Q 

charge  density  w.r.t.  cross-sectional  area 
(Cm-2) 

Q 

volumetric  charge  density  w.r.t.  total  volume 
(Cm-3) 

r 

reaction  rate  (mol  m-3  s-1) 

R 

ideal  gas  constant,  R  —  8.314  Jmol-1  K-1 
(J  mol-1  K-1) 

t 

time  (s) 

T 

temperature  (K) 

U 

voltage  (V) 

V 

velocity  (ms-1) 

V 

volume  (m3) 

V 

molar  volume  (m3  mol-1) 

(E 

)  atomic  diffusion  volumes  (cm3mol-1) 

w 

mass  fraction  (-) 

X 

mole  fraction  in  liquid-phase  (-) 

y 

mole  fraction  in  gas-phase  (-) 

z 

cell  coordinate  perpendicular  to  cell  plane  (m) 

z* 

number  of  transferred  electrons/single  charges 
(-) 

Greek 

symbols 

a 

heat  transfer  coefficient  (W  m-2  K-1) 

aa,  otc 

charge  transfer  coefficients  (anodic,  cathodic) 

(-) 

U'Vignes 

thermodynamic  factor  (Vignes  method)  (-) 

a' 

viscous  selectivity  factor  (-) 

K 

ratio  of  specific  heat  capacities,  chapter  A. 8  (-) 

r 

volumetric  charge  production  (Cm-3  s-1) 

e 

volume  fraction  (pore  volume  fraction  =  poros¬ 
ity)  (-) 

n 

overpotential  (V) 

rfn 

dynamic  viscosity  (Pa  s) 

X 

thermal  conductivity  coefficient  (W  m-1  K-1) 

A 

relative  water  content  in  membrane  (-) 

M 

chemical  potential  (J  mol-1) 

/xvis 

kinematic  viscosity  (mm2  s-1) 

V 

stoichiometric  coefficient  (-) 

p 

mass  density  (kg  m-3) 

r 

tortuosity  factor  (-) 

<P 

electrical  potential  (V) 

X 

non-ideality  coefficient  in  Flory-Huggins  ac¬ 
tivity  model  (-) 

Superscripts 

A 

anode  compartment  (supply  channel  structure) 

AC 

anode  catalyst  layer 

ACP 

polymer  phase  within  (AC) 

AD 

anode  diffusion  layer 

AF 

anode  feed 
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C  cathode  compartment  (supply  channel  struc¬ 

ture) 

CC  cathode  catalyst  layer 

CCP  polymer  phase  within  (CC) 

CD  cathode  diffusion  layer 

CF  cathode  feed 

eff  effective 

M  membrane  (PEM) 

vis  viscosity 

9  at  standard  conditions:  T6  —  298  K, 

pe  =  105  Pa 


Subscripts 

a 

anode 

air 

air 

BET 

BET  surface 

c 

cathode 

carbon 

carbon  material 

cat 

catalyst 

cat.layer  catalyst  layer 

cell 

cell 

cross 

crossover 

cu 

polymer  chain  unit 

CH3OH 

methanol 

C02 

carbon  dioxide 

dry 

dry 

eff 

effective 

eq 

equilibrium 

F 

feed 

gas 

gas 

(g) 

in  gas  state 

graphite 

graphite  material 

H+ 

proton 

h2o 

water 

i 

counting  index 

j 

counting  index 

Joule 

Joule  heating 

liquid 

liquid 

(1) 

in  liquid  state 

M 

solid  matrix  in  porous  materials 

n2 

nitrogen 

o2 

oxygen 

p 

at  constant  pressure 

pores 

in  pore(s) 

P 

polymer 

PTFE 

PTFE  (polytetrafluorethylene,  TEFLON™) 

PTFE-treated  Toray  TORAY™  carbon  paper,  treated 
with  PTFE 

sat 

saturated 

sound 

sound 

untreated  Toray  TORAY™  carbon  paper,  as  supplied 
by  manufacturer 

wet 

wet 

ter  could  e.g.  enable  reliable  dynamic  operation  with  a  min¬ 
imised  buffer  for  electrical  energy,  as  long  as  no  significantly 
improved  materials  for  the  DMFC  are  available.  Also,  clever 
control  strategies  could  possibly  enhance  the  power  density 
and  thus  the  fuel  efficiency  [1]. 

Both  tasks,  material  development  as  well  as  controller 
design,  can  only  be  effectively  addressed  if  realistic  math¬ 
ematical  process  models  are  available.  In  the  field  of  fuel 
cells  in  general,  a  number  of  publications  address  this  topic, 
but  mainly  in  the  area  of  hydrogen-fed  polymer  electrolyte 
membrane  fuel  cells  (PEMFC).  For  the  liquid-fed  DMFC  (1- 
DMFC)  only  a  small  number  of  mathematical  models  has 
been  published.  Table  1  shows  a  systematic  comparison  of 
recent  publications  (since  1997)  with  respect  to  the  key  fea¬ 
tures  of  mathematical  fuel  cell  models.  In  the  first  column  of 
the  table  the  respective  reference  numbers  are  given. 

The  second  column  of  Table  1  shows  whether  the  pre¬ 
sented  models  are  dynamic  (dyn.)  or  steady-state  (s.s.).  Here 
it  becomes  evident,  that  except  for  publications  from  our 
group  [1,2]  (also  including  this  paper),  all  published  mod¬ 
els  are  steady-state,  i.e.  they  are  used  to  predict  steady-state 
current  voltage  characteristics  and  concentration  profiles. 

The  next  block  of  seven  columns  presents  the  dimension¬ 
ality  of  the  models  in  the  seven  functional  layers  of  a  DMFC 
(A  =  anode  compartment/flow  channels,  AD  =  anode  diffu¬ 
sion  layer,  AC  =  anode  catalyst  layer,  M  =  membrane/PEM, 
CC  =  cathode  catalyst  layer,  CD  =  cathode  diffusion  layer, 
C  =  cathode  compartment/flow  channels).  A  blank  indicates 
that  this  part  of  the  DMFC  is  not  included  in  the  model,  a  zero 
means  that  this  element  is  described  by  a  lumped  parameter 
model  (usually  as  an  ideally  mixed-phase)  and  a  “1”  means 
that  this  element  is  modeled  one-dimensional.  In  case  of  the 
flow  channels  (A,C)  this  usually  means  along  the  channel, 
while  for  the  inner  layers  of  the  DMFC  this  means  perpen¬ 
dicular  to  the  cell  plane.  One  can  see  from  Table  1,  that  with 
respect  to  the  spatial  model  structure  the  different  models 
vary  significantly,  depending  on  the  focus  of  the  respective 
work.  A  number  of  models  does  not  account  for  the  cathode 
side  gas  transport  to  the  cathode  catalyst  layer,  assuming  this 
contribution  not  to  be  dominating  for  the  respective  operating 
conditions  [1, 3-7,2].  The  model  to  be  presented  in  this  work 
covers  all  structural  layers  of  the  DMFC. 

In  the  next  two  columns  of  Table  1 ,  the  type  of  the  applied 
electrode  kinetic  expressions  is  presented.  In  most  1-DMFC 
models  simple  Tafel  type  rate  expressions  are  applied,  only 
few  papers  use  Butler- Volmer  (BV)  type  expressions  [1,5,2] 
which  are  also  able  to  predict  open-circuit  overpotentials. 
Some  models  even  use  more  realistic,  complex  multi-step 
reaction  kinetics  for  the  electrochemical  methanol  oxidation 
[1,2,6,8-11], 

The  next  two  columns  of  Table  1  compare  the  descrip¬ 
tion  of  the  phase  situation  in  the  anode  and  cathode  pore 
structures.  It  is  well  known  that  in  the  1-DMFC  under  a  va¬ 
riety  of  operating  conditions,  on  the  anode  as  well  as  on  the 
cathode  side  a  liquid  (1)  and  a  gas  (g)  phase  can  coexist.  On 
the  anode  side  this  is  due  to  the  production  of  carbon  diox- 


Table  1 

Comparison  of  1-DMFC  models  published  since  1997 

Reference  dyn./s.s.  Dimensionality  Electrode  kinetics  Phase  situation  Mass  transport  models  Membrane  Energy 

swelling  balance 

(A)  (AD)  (AC)  (M)  (CC)  (CD)  (C)  Anode  Cathode  Anode  Cathode  Anode  Membrane  Cathode 
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ide  in  the  anode  catalyst  layer  and  the  low  solubility  of  this 
gas  in  liquid  water  methanol  solutions,  especially  at  elevated 
temperatures.  Therefore,  carbon  dioxide  bubbles  are  formed. 
Whether  this  takes  place  within  the  porous  catalyst  and  dif¬ 
fusion  layer  structures  or  only  within  the  flow  channels  is 
still  under  discussion,  but  a  variety  of  models  assumes  two- 
phase  flow  inside  the  anode  pore  structure  [4,5,8-13].  Inside 
the  cathode  pore  structure,  product  water  may  condense  and 
block  the  way  for  fresh  oxygen.  This  phenomenon  is  usually 
referred  to  as  cathode  flooding.  Also  here  it  is  not  fully  clear, 
whether  such  condensation  can  occur  within  the  (usually  hy¬ 
drophobic)  cathode  diffusion  layer,  or  only  on  the  surface 
of  the  diffusion  layer  inside  the  flow  channels.  Nonetheless, 
such  two-phase  behaviour  on  the  cathode  side  is  covered  by  a 
few  models  [8,12,13].  All  other  models  assume  pure  liquid- 
phase  on  the  anode  side,  and  pure  gas-phase  on  the  cathode 
side  of  the  1-DMFC. 

A  very  important  feature  of  each  1-DMFC  model  are  the 
chosen  mass  transport  descriptions  in  the  anode  and  cathode 
structures  and  inside  the  polymer  electrolyte  membrane.  Sev¬ 
eral  types  of  mass  transport  models  are  applied.  Simple  Fick 
diffusion  models  (F)  and  effective  Fick  models  (eff.F)  (us¬ 
ing  -  usually  experimentally  determined  -  effective  transport 
coefficients  instead  of  Fick  diffusivities)  do  not  account  for 
convective  flow  contributions  [1,2,4,13,14].  Therefore,  many 
models  feature  Nernst-Planck  (NP)  mass  transport  expres¬ 
sions,  which  combine  Fick  diffusion  with  a  superposed  con¬ 
vective  flow  [1-12,14].  The  latter  is  usually  calculated  from 
Darcy’s  law  using  different  formulations  of  the  hydraulic 
permeability  coefficient.  Instead  of  Darcy’s  law,  also  some 
models  use  Schlogl’s  formulation  (S)  for  the  convective  flow 

[1.2.5.12] .  This  also  accounts  for  electro-osmotic  flow  and 
can  thus  also  be  used  for  mass  transport  inside  the  PEM.  An 
alternative,  very  simple  way  of  incorporating  electro-osmotic 
flow  in  the  membrane  mass  transport  is  applying  so-called 
drag  coefficient  models  (dc)  which  assume  a  proportionality 
of  the  convective  water  and  methanol  flow  to  the  proton  flow 
[3,4,6,7,14].  The  last  popular  type  of  mass  transport  descrip¬ 
tion  is  the  Maxwell-Stefan  formulation  for  multi-component 
mixtures.  But  it  is  often  only  applied  to  gas-phase  transport 

[8.12] .  Only  one  model  so  far  (except  for  this  paper)  uses 
this  formulation  also  for  (liquid-phase)  mass  transport  in¬ 
side  the  PEM  [9-11].  Rarely  used  for  liquid  flow  are  surface 
diffusion  models  (50),  or  models  derived  from  irreversible 
thermodynamics  (irrev.TD)  [8].  All  mass  transport  models 
applying  effective  transport  coefficients  and  drag  coefficients 
(F,  eff.F,  NP,  S)  usually  only  yield  good  approximations  to 
experimental  data  for  a  very  limited  range  of  operating  con¬ 
ditions,  unless  the  coefficients  are  formulated  as  functions  of 
the  operating  conditions  (most  important  is  the  temperature 
dependence  of  all  mass  transport  parameters). 

The  second  last  column  of  Table  1  presents,  whether  the 
1-DMFC  models  account  for  varying  water  contents  inside 
the  PEM,  i.e.  the  swelling  behaviour  of  these  materials.  Most 
models  assume  a  fully  hydrated  PEM,  as  on  the  anode  side 
liquid  water  as  excess  component  is  present.  Only  few  mod- 
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els  skip  this  assumption.  In  one  case  [8]  the  water  uptake  of 
the  PEM  is  described  by  an  empiric  correlation  (developed 
from  experimental  data),  in  another  [9-11]  a  thorough  ther¬ 
modynamic  (TD)  model  is  formulated  based  on  change  of 
free  Gibbs  energy  inside  the  PEM  according  to  local  water 
content.  Also  in  this  paper,  a  thermodynamic  model  will  be 
presented  to  account  for  local  PEM  swelling. 

Finally,  the  last  column  of  Table  1  compares  the  published 
1-DMFC  models  with  respect  to  inclusion  of  energy  balances. 
Obviously  all  published  models  assume  an  isothermal  cell 
operation,  therefore  no  energy  balances  are  formulated.  The 
model  to  be  presented  in  this  paper  is  the  first  comprising  a 
full  energy  balance. 

Summing  up,  in  this  work  a  one -dimensional  model  of 
a  1-DMFC  will  be  presented,  which  is  different  from  so-far 
published  models  in  several  respects: 

•  The  model  is  dynamic,  allowing  also  to  predict  dynamic 
operation  (not  presented  in  this  paper). 

•  The  model  consequently  uses  a  Maxwell-Stefan  model 
for  all  types  of  mass  transport  in  all  functional  layers.  This 
enables  to  predict  mass  transport  correctly  for  a  vast  range 
of  operating  conditions  (especially  cell  temperatures  be¬ 
tween  ambient  and  90  °C)  with  one  single  set  of  mass 
transport  parameters. 

•  Membrane  water  uptake  and  (local)  swelling  behaviour 
are  accounted  for  by  applying  a  Flory-Huggins  activity 
model. 

•  Heat  conduction  and  local  heat  production  are  fully  ac¬ 
counted  for  by  a  complete  energy  balancing. 

The  mean  simplifications  of  the  model  are: 

•  Assumption  of  pure  liquid-phase  in  the  anode  structures 
and  pure  gas  phase  in  the  cathode  structures. 

•  No  spatial  discretisation  of  catalyst  layers. 

•  Application  of  Butler- Volmer  type  rate  expressions  for 
both  electrode  reactions. 

To  evaluate  the  model,  experiments  were  performed  in 
which  methanol  and  water  crossover  fluxes  through  the  PEM 
were  measured  (together  with  the  cell  voltage)  under  a  variety 
of  operating  conditions  (anode  feed  temperatures  from  30  to 
90  °C,  full  range  of  cell  current  densities). 

In  this  paper,  the  focus  will  be  on  describing  the  struc¬ 
ture  of  the  model  and  presenting  the  results.  More  details, 
especially  with  respect  to  the  mass  transport  and  the  Flory- 
Huggins  activity  model  of  the  PEM,  phase  equilibria  between 
PEM  and  the  catalyst  layers,  the  derivation  of  all  model  pa¬ 
rameters,  as  well  as  a  number  of  additional  experimental  re¬ 
sults  can  be  found  in  [15]. 

2.  Experiments 

To  be  able  to  measure  methanol  and  water  crossover  fluxes 
in  a  DMFC,  a  fully  automated  miniplant  was  constructed  and 
an  own  DMFC  design  was  developed  in  order  to  be  able  to 


influence  all  materials  and  structures  within  the  DMFC.  In 
the  following,  the  miniplant  and  the  in-house  DMFC  design 
will  be  shortly  presented. 

2.1.  Applied  in-house  DMFC  design 

The  experiments  were  carried  out  using  a  single  cell 
DMFC  fed  with  air  and  liquid-methanol-water  solutions.  A 
detailed  description  of  the  DMFC  design  can  be  found  in 

[15]. 

The  identical  anode  and  cathode  monopolar  plates  are 
made  from  graphite  material  (thickness  7  mm,  material  code 
FU4369)  supplied  by  SchunkKohlenstofftechnik  (Germany). 
The  necessary  flowbed  structures  for  the  reactant  distribution 
over  the  MEA  surface  are  millcut  into  the  plates  (Fig.  1).  They 
consist  of  parallel  channels  of  2  mm  width  and  2  mm  depth, 
with  1mm  wide  ribs  between  them.  A  distributor  (Fig.  1, 
top)  and  collector  channel  (Fig.  1,  bottom)  connect  the  par¬ 
allel  channels  to  the  inlet  and  outlet  ports,  respectively.  The 
media  (air  and  methanol-water  solution)  are  supplied  in  one 
corner  of  the  rectangular  flowbed  and  leave  at  the  opposite 
corner  (flow  direction  indicated  in  Fig.  1).  The  flowbed  itself 
has  the  outer  dimensions  65  mm  x  40  mm,  identical  to  the 
catalyst  layer  on  the  MEAs,  which  leads  to  an  active  area  of 
As  =  26  cm2. 

As  diffusion  layers  PTFE-coated  TORAY  carbon  paper 
(TGP-H-060)  is  used,  with  a  PTFE  loading  between  20  and 
25  mass%  with  respect  to  the  uncoated  material. 

Finally,  the  membrane  electrode  assemblies  (MEA)  are 
prepared  fromNAFION™  N-105  membrane  foil,  onto  which 
the  catalyst  layers  are  applied  using  an  airbrush  technique 
refined  by  ZSW  Ulm  (Germany)  [16].  The  anode  catalyst 
layer  features  a  catalyst  loading  of  5  mg  cm-2  (unsupported) 
platinum  ruthenium  black  (Alfa  Aesar  Johnson  Matthey 
HiSPEC™  6000)  and  a  NAFION™  content  of  15mass% 
relative  to  the  metal  loading  (i.e.  0.75  mg  cm-2).  The  cathode 
catalyst  layer  has  the  same  metal  loading,  but  as  catalyst 


Fig.  1 .  Photo  of  monopolar  plate  showing  inlet,  outlet  and  flow  direction  in 
flowbed  structure. 
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(unsupported)  platinum  black  is  used  (Alfa  Aesar  Johnson 
Matthey  HiSPEC™  1000)  and  the  NAFION™  content  is 
10mass%  relative  to  the  metal  loading  (i.e.  0.5  mg  cm-2). 

The  DMFC  is  completed  by  gold-plated  copper  plates  as 
current  collectors  and  stainless  steel  plates  for  bracing  the 
whole  sandwich  structure.  A  torque  of  5  Nm  is  exerted  on 
the  screws,  which  hold  together  the  steel  back  plates.  After 
assembly,  each  DMFC  is  conditioned  and  evaluated  by  oper¬ 
ation  with  pure  humidified  hydrogen  and  air  for  three  times 
8  h,  before  it  is  operated  on  methanol  solutions. 

2.2.  Experimental  setup 

For  automated  testing  of  DMFCs,  a  miniplant  was  de¬ 
signed  using  the  process  control  system  PC-S7/WinCC  by 
Siemens.  It  enables  automatic  testing  procedures,  with  a  spe¬ 
cial  focus  on  dynamic  operation.  Fig.  2  shows  a  simplified 
flowsheet  of  the  miniplant.  All  details  can  be  found  in  [15]. 

The  DMFC  cathode  is  supplied  with  dry  air  (dew  point 
~3  °C)  at  flow  rates  between  0.4  and  5.0scbmh-1  (mass 
flow  controller  F101,  type  Mass6020  by  Burkert  AG,  Ger¬ 
many)  at  cathode  outlet  pressures  of  ambient  up  to  5  bars 
absolute  (1.5  x  105  Pa).  The  air  is  pre-heated  in  a  plate  heat 
exchanger  (W101),  air  temperatures  and  pressures  are  mea¬ 
sured  at  the  cathode  inlet  and  outlet.  At  the  cathode  outlet, 


also  the  relative  humidity  of  the  air  is  measured  (Q202,  type 
HygroClip  IE  by  rotronic  AG,  Switzerland).  Finally  the  cath¬ 
ode  exhaust  air  enters  a  condensor,  where  it  is  dried  to  reach  a 
dew  point  below  10  °C  (condensate  is  collected).  The  dry  air 
is  sent  into  a  fume  hood,  while  its  oxygen  and  carbon  diox¬ 
ide  contents  are  measured.  The  oxygen  sensor  (Q204)  is  a 
paramagnetic  sensor  (PAROX  1000  H  by  MBE  AG,  Switzer¬ 
land),  while  carbon  dioxide  is  measured  using  an  FT-IR- 
sensor  (Q203,  type  OEM-NDIR  EGC-5%  by  Pewatron  AG, 
Switzerland). 

On  the  DMFC  anode  side,  a  liquid  recycle  loop  is  installed. 
It  consists  of  two  alternative  cycles,  one  for  methanol-water 
solution  and  one  for  pure  water.  The  purpose  of  this  is  to 
enable  a  stepped  or  pulsed  periodic  operation  of  the  DMFC, 
where  the  anode  feed  is  changed  stepwise  between  methanol- 
water  solution  and  pure  water  automatically.  Both  branches  of 
the  anode  cycle  feature  vessels  for  pressure  equilibration  and 
carbon  dioxide  removal  (B 1  and  B2),  gear  pumps  (P401  and 
P402)  and  heat  exchangers  (W403  and  W406).  Flow  rates 
between  0.3  and  5  dm3  min-1  can  be  achieved.  Automatic 
valves  (V403/V404  and  V408/V409)  enable  a  flexible  and 
practically  immediate  change  between  methanol-water  so¬ 
lution  and  pure  water  anode  feed  without  causing  significant 
disturbances  in  liquid  flow  rate  and  pressure.  The  flow  rate 
is  measured  by  a  Coriolis-type  sensor  (F401,  type  MASS 


Fig.  2.  Flow  sheet  (simplified)  of  DMFC  miniplant  with  major  components  and  sensors. 
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2100  DI6  by  Danfoss,  Denmark).  Like  on  the  cathode  side, 
the  medium  temperature  and  pressure  are  measured  at  the 
cell  inlet  and  outlet.  The  methanol  concentration  of  the  an¬ 
ode  inlet  medium  is  measured  online  using  an  ultrasound 
sensor  (Q401,  type  LiquiSonic30  by  SensoTech  GmbH,  Ger¬ 
many),  based  on  the  influence  of  the  methanol  concentra¬ 
tion  on  the  speed  of  sound  in  methanol-water  solutions. 
This  sensor  is  used  in  a  methanol  concentration  controller, 
which  as  actuators  uses  two  dosing  pumps  for  pure  methanol 
(P351)  and  pure  water  (P301)  (mzr-2905  by  HNP  Mikrosys- 
teme,  Germany).  Methanol  concentrations  between  0  and 
1.5  mol  dm-3  can  be  detected  and  controlled.  The  flow  rates 
of  the  dosing  pumps  can  be  controlled  in  the  range  from  0.2 
up  to  18  cm3  min-1.  To  adjust  the  anode  pressure  and  also 
to  strip  off  carbon  dioxide,  the  recycle  vessels  are  equipped 
with  a  nitrogen  purge/blanket.  The  anode  pressure  can  be 
controlled  in  the  range  between  ambient  and  5  bars  abso¬ 
lute  (1.5  x  105  Pa).  The  liquid  inlet  temperature  (which  is 
also  the  DMFC  temperature  due  to  the  applied  high  flow 
rates)  can  be  controlled  in  the  range  between  —20  and  +90  °C 
(253,  ...,363 K). 

The  DMFC  is  electrically  connected  to  a  potentiostat 
(HP60-50  by  Wenking  GmbH,  Germany),  which  enables  op¬ 
eration  of  fuel  cells  from  below  1  W  up  to  1  kW  at  a  maximum 
current  of  50  A.  Galvanostatic  as  well  as  potentiostatic  op¬ 
eration  is  possible,  with  the  possibility  to  automatically  run 
user-defined  load  scenarios. 

The  array  of  sensors  around  the  DMFC,  i.e.  potentiostat, 
flow  rates,  inlet  and  outlet  concentrations  of  key  components 
(oxygen,  methanol,  water,  carbon  dioxide)  enables  full  online 
material  balancing  of  these  components.  Assuming  full  and 
immediate  oxidation  of  crossover  methanol  on  the  cathode 
catalyst,  from  the  sensor  information  also  the  methanol  and 
water  crossover  fluxes  from  anode  to  cathode  (i.e.  through 
the  PEM)  can  be  calculated  (for  details  see  [15]). 


2.3.  Performed  experiments 

The  online  balancing  function  of  the  miniplant  allows  to 
measure  the  steady-state  methanol  and  water  crossover  fluxes 
through  the  PEM  of  the  DMFC.  These  crossover  fluxes  (to¬ 
gether  with  the  current-voltage  characteristics)  were  mea¬ 
sured  for  a  broad  range  of  cell  temperatures  (30-90  °C)  as 
functions  of  the  cell  current  density.  The  experimental  results 
are  presented  in  chapter  4  (model  validation)  in  comparison 
to  simulation  data.  All  operating  conditions  are  given  there. 


3.  Model  formulation 

In  this  chapter,  a  dynamic  model  of  a  single  cell  DMFC 
will  be  presented.  This  model  represents  the  cross-sectional 
structure  of  the  DMFC,  which  is  depicted  in  Fig.  3.  The  model 
is  one-dimensional,  perpendicular  to  the  cross-sectional  area 
of  the  cell.  All  state  variables  are  assumed  to  be  constant  in  the 
other  two  space  coordinates.  As  can  be  seen  in  Fig.  3,  on  this 
level  of  decomposition  the  cell  consists  of  seven  sequentially 
connected  phases:  anode  and  cathode  compartments  (i.e.  the 
flowbeds),  both  diffusion  and  catalyst  layers  and  the  PEM  in 
the  middle. 

3.1.  Basic  model  assumptions 

The  basic  model  assumptions,  which  apply  to  all  elements, 
are: 

•  all  gas  phases  obey  the  ideal  gas  law, 

•  all  inert  gases  (components  of  air,  like  nitrogen,  argon  etc.) 
are  merged  to  “nitrogen”, 

•  in  liquid  phases  nonideal  mixing  behaviour  is  not 
accounted  for,  variations  of  activity  coefficients  are 


process  unit  level: 


phase  level: 

anode 

compart- 

i  "" 


anode  . 
diffusion 


anode 

catalyst 


proton 

exchange; 


cathode 

catalyst 


cathode 

diffusion 


cathode 
compart- 
'  ± 


phase  boundaries 
(internal  boundary  conditions) 


Fig.  3.  Schematic  of  DMFC  layers,  component  mass  fluxes  and  model  structure. 
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neglected  due  to  the  strong  dilution  of  methanol  in  water 
(this  assumption  is  not  valid  in  the  PEM), 

•  in  each  element,  only  the  existence  of  one  thermodynamic 
phase  is  assumed,  i.e.  on  the  anode  side  formation  of 
carbon  dioxide  bubbles  is  neglected,  on  the  cathode 
side  the  formation  of  a  liquid-phase  due  to  condensation 
(so-called  cathode  flooding), 

•  all  charge  balances  are  quasi-stationary, 

•  electroneutrality  is  assumed  in  the  PEM, 

•  and  finally,  ohmic  losses  other  than  in  the  PEM  are 
neglected. 

To  assume  ideal  gas  behaviour  seems  justified,  as  the  max¬ 
imum  system  pressures  do  not  exceed  5  bars  (5  x  105Pa) 
and  the  temperature  range  does  not  exceed  10,  ... ,  90  °C 
(283, . . . ,  363  K).  Also  all  components  (methanol,  water, 
oxygen,  carbon  dioxide  and  nitrogen)  have  a  low  molecu¬ 
lar  weight. 

The  assumption  of  solely  single-phase  behaviour  is  the 
major  simplification  of  the  model.  Under  a  variety  of  prac¬ 
tical  operating  conditions  it  has  been  observed,  that  liquid 
water  forms  in  the  cathode  compartment  and  gas  bubbles  of 
carbon  dioxide  form  in  the  anode  compartment.  The  question 
which  phase  situation  is  present  within  the  porous  diffusion 
and  catalyst  layers,  though,  has  not  yet  been  answered  satis¬ 
factorily. 

The  final  assumption  that  no  ohmic  drops  are  accounted 
for  other  than  that  in  the  PEM  is  based  on  the  fact  that  all  elec¬ 
tron  conducting  parts  of  the  DMFC  (bipolar  plates,  diffusion 
layers,  catalyst  layers)  are  made  from  very  good  electron  con¬ 
ducting  materials  (graphite,  noble  metals),  and  that  especially 
the  diffusion  and  catalyst  layers  are  very  thin,  while  having 
a  large  contact  area. 


3.2.  Mass  transport  and  balancing 


For  the  description  of  mass  transport  in  porous  structures, 
the  generalised  Maxwell-Stefan  approach  is  used,  in  the  for¬ 
mulation  proposed  by  Krishna  and  Wesselingh  [17,18].  It  is 
based  on  a  mechanical  equilibrium  between  driving  forces 
acting  on  a  species  j  and  friction  forces  between  this  species 
and  all  other  species  i  around  it: 
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On  the  left  hand  side  of  Eq.  (1)  are  four  terms  describing 
the  driving  forces.  Term  one  and  two  describe  diffusive  driv¬ 
ing  forces  resulting  from  gradients  in  the  chemical  potentials 
p.j  (Jmol-1)  (first  term:  at  constant  pressure,  second  term: 
pressure  influence),  term  three  is  the  driving  force  due  to  a 
superficial  viscous  flow  resulting  from  a  gradient  in  the  total 
pressure p  (Pa)  and  the  fourth  term  represents  the  electrostatic 
force  acting  on  charged  species  due  to  a  gradient  of  the  elec¬ 


tric  potential  field  cp  (V)  (migration).  On  the  right  hand  side 
are  two  terms  describing  friction  forces:  the  sum  accounts  for 
the  friction  between  species  j  and  all  other  mobile  species  i 
(x  are  mole  fractions,  n  are  flux  densities  in  (mol  m-2  s-1)), 
and  the  second  term  represents  friction  between  species  j  and 
the  (stationary)  solid  matrix  (lower  index  “M”). 

The  most  important  parameters  in  this  equation  are  the  bi¬ 
nary  Maxwell-Stefan  diffusion  coefficients  D  (m2  s-1).  The 
lower  indices  denote  the  two  respective  species,  an  upper  in¬ 
dex  “eff”  means  that  this  is  an  effective  diffusion  coefficient 
taking  into  account  the  porosity  s  (-)  and  the  tortuosity  r  (-) 
of  the  solid  matrix,  while  those  diffusion  coefficients  without 
the  index  “eff”  are  valid  for  free  space  binary  interactions. 
For  a  more  detailed  treatment  of  the  binary  diffusion  coef¬ 
ficients  refer  to  Appendix  A. 9.  The  other  parameters  of  the 
presented  form  of  the  generalised  Maxwell-Stefan  equations 
are  explained  in  the  respective  following  sections  and  in  the 
list  of  symbols. 

The  modeling  concept  is  based  upon  a  finite  volume 
discretisation  along  only  one  spatial  coordinate  z.  For  finite 
volume  element  simulations  it  is  well  known,  that  the  simul¬ 
taneous  treatment  of  diffusive  fluxes  and  convective  flow  can 
lead  to  numerical  problems  if  the  viscous  flow  contribution 
in  relation  to  the  overall  mass  transport  is  high.  To  prevent 
such  problems,  the  viscous  flow  term  (term  three  on  the  left 
hand  side)  of  Eq.  (1)  is  skipped  and  the  Maxwell-Stefan 
equation  is  formulated  only  for  the  individual  driving  forces 
acting  on  the  species  j  (diffusion  and  migration) 
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In  Eq.  (2)  in  the  friction  terms  on  the  right  hand  side  not 
the  overall  flux  densities  n  j  appear  but  only  the  individual 
flux  densities  j  j. 

The  overall  flux  densities  are  then  calculated  from  the 
individual  flux  densities  and  an  additive  term  for  the  pressure- 
driven  convective  contribution: 


"/  =  Jj  +  CjVp. 


(3) 


Here  not  the  molar  concentration  with  respect  to  the  vol¬ 
ume  of  the  fluid  phase,  c  (mol  m-3),  is  used,  but  a  modified 
concentration  c  (molm-3)  with  respect  to  the  total  volume 
including  the  porous  matrix.  i;p  is  the  convective  velocity 
(ms-1),  which  is  a  function  of  the  total  pressure  gradient 
according  to  Darcy’s  law: 
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The  parameters  are  the  dynamic  viscosity  if'1'  (Pa  s)  of  the 
mixture  and  the  hydraulic  permeability  coefficient  B°  (m2). 
The  latter  has  to  be  determined  experimentally,  but  for  some 
simple  geometries  correlations  are  known.  The  simplest  case 
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is  the  flow  through  parallel  straight  tubular  pores  (Poiseuille 
flow): 


So  = 


32 


(5) 


In  the  following,  this  (simplifying)  approach  will  be  used. 

In  the  implementation  of  the  transport  equations  for  a  spa¬ 
tially  discretised  model,  an  upwind  scheme  is  used  which 
always  uses  the  “upwind”  concentrations,  i.e.  the  concen¬ 
trations  in  the  left  or  the  right  neighbouring  control  vol¬ 
ume  depending  on  the  direction  of  the  convective  flux, 
to  calculate  the  convective  contribution  to  the  molar  flux 
densities. 

Using  the  overall  flux  densities,  the  general  form  of  the 
component  mass  balances  is 
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with  the  reaction  rates  r \  of  all  occurring  (electro-)  chemical 
reactions  (molm-3  s-1)  and  the  stoichiometric  coefficients 
vj,k  (-)  of  component  i  in  reaction  k. 

Finally  a  total  mass  balance  can  be  formulated  based  on 
the  continuity  equation: 


3.3.  Energy  transport  and  balancing 

Within  the  DMFC,  not  only  a  variety  of  different  mass 
transport  phenomena  occur  simultaneously.  The  same  is  true 
for  energy  transport  and  production.  As  the  DMFC  consists 
of  porous  layers  in  which  mobile  species  are  transported,  en¬ 
ergy  transport  can  take  place  both  due  to  transport  bound  to 
the  moving  species  and  thermal  conduction.  The  latter  takes 
place  in  the  mobile  phase  as  well  as  in  the  stationary  solid 
matrices  in  the  different  layers.  Additionally,  chemical  and 
electrochemical  reactions  take  place  in  the  catalyst  layers 
(AC)  and  (CC),  and  finally  within  the  membrane  (M)  a  spa¬ 
tially  distributed  heat  production,  Joule  heating,  occurs,  due 
to  the  transport  of  charged  species  in  an  electric  field. 

To  get  a  most  simple  model  description  of  all  these  phe¬ 
nomena,  it  makes  sense  to  formulate  two  independent  energy 
flux  densities:  enthalpy  flux  densities 

e  =  5>;/U(7’)  (12) 

j  j 

(Jm“2  s-1),  which  are  coupled  to  the  mass  flux  densities 
n;  (molm-2s-1)  of  the  mobile  species  and  their  specific 
enthalpies  h;  (J  mol-1),  and  heat  flux  densities 
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where  p  (kgm-3)  is  the  fluid  mixture  density  and  v  (ms-1) 
is  the  mean  (superficial)  mixture  velocity.  For  compressible 
fluids  the  relation  between  pressure  and  density  at  constant 
entropy  (which  is  the  case  for  moderate  pressures  as  typical 
for  DMFC  operation)  is  given  by 

(!)s="»-  ,8» 

where  i>SOund  (ms-1)  is  the  speed  of  sound  in  the  fluid  (see 
Appendix  A)  and  p  (Pa)  is  the  local  pressure. 

Finally,  the  product  of  mean  mixture  density  and  mean 
velocity  is  the  total  mass  flux  density  (kgm-2  s-1): 


(J  m-2  s-1)  due  to  thermal  conduction  (Fourier  law). 

In  Eq.  (13)  /.dl  (W  m-1  K-1)  stands  for  the  local  effec¬ 
tive  thermal  conductivity  coefficient.  The  upper  index  “eff” 
denotes  that  it  is  dealt  with  a  mixture  of  a  fluid  and  a  solid- 
phase,  which  both  contribute  to  the  thermal  conduction.  The 
effective  thermal  conductivity  has  to  be  calculated  from  the 
thermal  conductivities  of  both  phases  taking  into  account 
their  volume  fractions.  These  calculations  can  be  found  in 
Appendix  C. 

The  specific  enthalpies  of  the  mobile  species  are  calcu¬ 
lated  from  the  specific  enthalpies  of  formation,  A -pH/,  and 
the  mean  heat  capacities  Cpj: 


m  tot  =  pv 


(9) 


Combining  Eqs.  (7) — (9)  one  can  formulate  the  total  mass 
balance  in  terms  of  the  pressure  as  variable,  and  its  time 
derivative  as  a  function  of  the  total  mass  flux  density: 


dp  _  2  9m  tot 

~dt  ~  ^sound  ’ 


(10) 


In  the  simulation,  the  total  mass  flux  density  is  cal¬ 
culated  simply  as  the  sum  of  the  component  mass  flux 
densities,  which  in  turn  are  the  products  of  the  molar 
flux  densities  and  the  respective  molecular  weights  Mj 
(kg  mol-1): 


m  tot  =  njMj- 

j 


(11) 


h  j  =  A ¥Hej  +  j  Cpj(T)  dT  «  AF H]  +  (T  -  Te)  •  CpJ 

(14) 

By  using  enthalpy  fluxes,  all  heats  of  reactions  are  ac¬ 
counted  for  automatically  without  the  need  for  a  heat  pro¬ 
duction  term  in  the  energy  balances  of  the  respective  DMFC 
layers. 

Finally,  Joule  heating  ejouie  (W  m-3)  due  to  charge  trans¬ 
port  is  described  by  the  general  equation 

.  dtp 

Ctoule  —  i  — ,  (15) 

dz 

(W  m-3).  It  is  always  positive  and  independent  of  the  direc¬ 
tion  of  the  charge  flux. 
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Combining  all  three  energy  flux  densities  yields  the  gen¬ 
eral  energy  balance: 


dT 

Ik 


1 

(pc  p) 


de 

3  z 


3  q 

„  h  eJoule 


(16) 


The  local  effective  volumetric  heat  capacities  ( 'pcp )  are 
caclulated  from  the  local  porosities,  concentrations  and  heat 
capacities  of  the  present  components.  The  calculation  is  pre¬ 
sented  in  Appendix  D. 


3.4.  Charge  transport  and  balancing 


In  the  PEM  no  charge  production  occurs  and  local  elec¬ 
troneutrality  is  assumed,  therefore  one  ends  up  with  the  quasi¬ 
stationary  charge  balance: 


3  z 


(22) 


3.5.  Overpotentials  and  cell  voltage 


For  the  formulation  of  rate  equations  for  the  electrochem¬ 
ical  electrode  reactions,  one  needs  a  definition  for  the  elec¬ 
trode  overpotentials,  which  is  given  by: 


Charge  transport  in  the  DMFC  is  bound  to  protons  within 
the  membrane  material  (proton  conductor)  and  electrons 
in  the  electrical  circuit  (electron  conductor).  Therefore,  the 
charge  flux  density  i  (Am-2)  within  the  proton  conductor  is 
coupled  to  the  molar  flux  density  of  protons  by  Faradays  law 
(z*  =  number  of  single  charges  exchanged  per  molecule): 

'  =  Z*H+FnH+  =  FnH+.  (17) 

Charge  production/consumption  due  to  electrochemical 
reactions  is  similarly  linked  to  the  reaction  rates  with  respect 
to  the  pore  volume,  r^  (molm-3  s-1),  of  the  respective  elec¬ 
trochemical  reactions  (as  function  of  the  respective  overpo¬ 
tential  r\  (V)): 

y  pores 

ik(ri)  =  Zn+FvH+,krk(.r})—^-.  (18) 

In  both,  the  membrane  and  the  catalyst  layers,  charge  bal¬ 
ances  have  to  be  formulated  accounting  for  the  outer  electri¬ 
cal  cell  current  density  iceii  (Am-2)  as  well  as  for  the  above 
mentioned  current  densities  resulting  from  proton  flux  in  the 
membrane  and  the  charge  production  in  both  catalyst  layers. 
The  general  form  of  the  charge  balances  is 


ri=  A0  -  A </>f=0  (23) 

The  overpotentials  ?;  (V)  are  defined  as  the  difference  be¬ 
tween  the  real  electrode  potential  A(p  (V)  (w.r.t.  standard  hy¬ 
drogen  electrode)  and  that  at  open  circuit  condition  (i.e.  no 
cell  current,  i  =  0)  and  thermodynamic  standard  conditions 
(pressure  105Pa,  temperature  298  K,  all  reactants  activities 
equal  one,  upper  index  0).  For  the  DMFC  the  following  values 
can  be  found  in  the  literature: 


Anodic  methanol  oxidation: 

A<i=0  =  0.02  V 

(24) 

Cathodic  oxygen  reduction: 

A  <i=0  =  1-23V 

(25) 

Then  the  cell  voltage  can  be  calculated  from  the  reversible 

open  circuit  cell  voltage  at  the  above  mentioned  standard  con¬ 
ditions  (  U&  n  /=0  ~  1.21V),  the  anode  and  cathode  overpo¬ 
tentials,  ?7a  and  t]c  (V),  respectively,  and  the  Ohmic  losses 
within  the  membrane  represented  by  the  total  difference  in 
the  polymer-phase  electrical  potentials  A <pM  (V): 

(/cell  =  (/cell, i=0  -  7a  +  7c  -  A </>M.  (26) 


a Q 

3 1 


3  i 
3  z 


+  n  7) 


(19) 


where  Q  is  the  volumetric  charge  density  (Cm-3)  and  /’ 
is  the  volumetric  charge  production  (C  m-3  s- 1 )  by  electro¬ 
chemical  reactions: 


H,)  =  $>(,).  (20) 

k 

Simplifyingly  it  is  assumed  that  the  charge  balances  are 
fast  compared  to  all  other  balances  (material  and  energy) 
which  leads  to  quasi-stationary  formulations.  In  the  catalyst 
layers,  when  balancing  the  electrons,  only  the  charge  flux  to 
or  from  the  adjacent  diffusion  layer  (electron  conductor)  oc¬ 
curs,  which  is  the  cell  current  density  ;cen .  Within  the  catalyst 
layers  electrons  are  produced  or  consumed  by  the  electro¬ 
chemical  reactions  as  described  above.  Therefore,  the  final 
quasi-stationary  charge  balances  have  the  form: 

0  =  (cell  -£>(7).  (21) 

k 


3.6.  Anode  compartment  (A) 


The  anode  compartment  is  assumed  to  be  a  spatially  con¬ 
centrated  phase  element  (ideally  mixed,  CSTR  behaviour). 
It  has  one  inlet  (feed,  index  AF)  and  one  outlet,  and  it  is 
connected  to  the  anode  diffusion  layer  (AD).  Fresh  water 
methanol  solution  is  fed  at  the  inlet,  with  a  supposedly  very 
low  carbon  dioxide  content.  Methanol  and  water  are  trans¬ 
ported  through  the  anode  diffusion  layer  towards  the  anode 
catalyst,  while  produced  carbon  dioxide  is  transported  in 
the  opposite  direction  to  leave  the  diffusion  layer  into  (A), 
and  on  this  way  being  removed  through  the  anode  outlet 
(see  Fig.  3).  The  material  balances  are: 


1 

pA 

1 

pA 


(TAFcf  -  Faca  +  As«ad(,ad  =  Q)) 

(Taf(caf  -  cf)  +  AstiAD(zAD  =  0)) 


(27) 


with  j  =  H20,  CH3OH,  C02.  In  Eq.  (27)  Faf  is  the  feed 
volume  flow  rate  (m3s-1),  PA  is  the  total  volume  of  the 
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channels  (m3),  cAF  are  the  feed  concentrations  (mol  m  3)  and 
FA  is  the  outlet  volume  flow  rate  (m3  s_1).  The  difference 
between  the  two  flow  rates  Faf  and  FA  can  be  assumed  to  be 
small,  therefore  it  is  neglected  which  results  in  the  simplified 
formulation  given  in  the  second  line  of  Eq.  (27). 

The  anode  pressure  pA  as  well  as  the  temperature  TA 
are  given  as  input  parameters  which  are  known  from  the  ex¬ 
periments.  Therefore,  no  total  mass  and  energy  balances  are 
formulated  here. 


3. 7.  Cathode  compartment  ( C) 


The  cathode  compartment’s  structure  is  similar  to  that  of 
the  anode  compartment.  The  inlet  is  fed  with  air  with  the  com¬ 
ponents  oxygen,  nitrogen  (including  all  other  inert  gases), 
water  and  carbon  dioxide.  Oxygen  enters  the  cathode  diffu¬ 
sion  layer  (CD)  while  water  and  carbon  dioxide  produced 
in  the  cathode  catalyst  layer  (CC)  are  transported  out  of  the 
DMFC  (see  Fig.  3) 

The  material  balances  are: 

^  =  ^(FCFcf  -  Fcf  +  AsnJ] D(zCD  =  dCD))  (28) 

with  j  =  N2,  O2,  H2O,  CO2. 

Here  Fcf  is  the  feed  volume  flow  rate  (m3  s_l),  Vc  is 
the  total  volume  of  the  channels  (m3),  cf  are  the  feed  con¬ 
centrations  (mol  m  3  )  and  Fc  is  the  outlet  volume  flow  rate 
(m3  s“ 1 ).  The  latter  can  be  calculated  from  the  inlet  flow  rate 
and  the  molar  fluxes  exchanged  with  the  diffusion  layer  (CD) 
(quasi-stationary  total  material  balance): 


pC  _  pCF 


RT 


c  »(zCD  =  dCD) 


(29) 


with  j  =  N2,  O2,  H2O,  CO2. 

The  cathode  pressure  pc  as  well  as  the  gas  and  the  bipolar 
plate  temperatures,  Tc  and  73arbon ,  respectively  (explanation 
for  two  temperatures  in  Section  3.9),  are  input  parameters 
known  from  experiments.  Therefore,  as  in  (A),  no  energy 
and  no  total  mass  balances  are  necessary. 


3.8.  Anode  diffusion  layer  (AD) 

The  anode  diffusion  layer  connects  the  anode  compart¬ 
ment  (A)  and  the  anode  catalyst  layer  (AC)  (see  Fig.  3).  It  con¬ 
sists  of  a  chemically  inert  carbon  fibre  material  coated  with 
a  certain  amount  of  PTFE  (20-25  wt.%).  It  supplies  educts 
(methanol  and  water)  to  the  anode  catalyst  and  removes  the 
carbon  dioxide  from  there.  It  also  collects  the  electrons  from 
the  anode  reaction  and  ensures  good  electric  contact  with 
the  bipolar  plate  (current  collector).  All  fluxes  are  assumed 
to  occur  only  perpendicular  to  the  surface  plane  (i.e.  in  z- 
direction).  Data  about  porosity  and  other  physical  parameters 
are  given  in  Appendix  B. 


It  is  assumed  that  the  complete  pore  space  is  filled  with 
liquid  methanol-water  solution  and  soluted  carbon  dioxide. 
Formation  of  a  gas  phase  from  carbon  dioxide  is  neglected 
(see  assumptions  in  Section  3.1). 

For  the  three  species  water  (H2O),  methanol  (CH3OH) 
and  carbon  dioxide  (CO2),  the  material  balances  (formulated 
in  molar  concentrations)  are: 


dcf 

dt 


1 


dnf 


c-AD 


dz 


(30) 


‘'pores 

with  j  =  H20,  CH3OH,  C02 

The  total  mass  balance,  analogue  to  Eq.  (10),  is  given  by 


dp 


,AD 


AD 


dt 


_/  AD  s2  9wtot 
'  sound/ 


Finally,  the  energy  balance  is 


dTAD 

dt 


1 

(ipf,)) 


,AD 


3eAD 

dz 


dqAD 

dz 


(31) 


(32) 


analogue  to  Eq.  (16),  but  Joule  heating  due  to  electron  trans¬ 
port  is  neglected  as  the  Ohmic  resistance  (and  therefore  the 
electric  potential  gradient)  in  the  carbon  paper  can  be  as¬ 
sumed  to  be  negligible. 


In  Eq.  (32)  the  molar  mixture  heat  capacity  (pc„)AD 


is  ap- 
AD 


proximatedby  the  value  for  pure  water  (see  Appendix  D).  e 
represents  the  sum  of  the  enthalpy  flux  densities  connected 
with  the  material  fluxes  and  qAD  is  the  heat  flux  density  due 
to  thermal  conduction  in  the  liquid  mixture  and  the  pore  walls 
(carbon  fibres). 

For  constant  activity  coefficients  and  pure  liquid-phase  the 
Maxwell-Stefan  mass  transport  equation,  Eq.  (2).  simplifies 
to: 


3vad 
„ad  axj 

tot  dz 


=  £ 


AD  :AD  _  AD  -AD 
■*/  Jj  xj  h 


■AD 


D 


AD.eff 

ij 


D 


AD  ,eff 
JM 


(33) 


with  j  =  H20,  CH3OH,  COi. 

~A] 
Hot 

species  concentrations: 


The  total  concentration  cA?  (mol  m  3)  is  the  sum  of  all 


AD 


J2cjd 


(34) 


with  j  =  H20,  CH3OH,  C02 

As  demonstrated  by  Krishna  and  Wesselingh  [18],  this 
flux-implicit  set  of  transport  equations  can  be  transformed  to 
get  explicit  formulations  for  the  flux  densities: 


0'AD)  = 


ADroADi-I, 


where  the  elements  of  the  transport  matrix  [Z?AD]  are: 


(35) 


diagonal  elements  :  BAD  — 


1 

r-.  AD.eff 
vi  M 


..AD 

EXk 

n AD.eff  ’ 
kjti  V ik 


(36) 


,.AD 


all  other  elements  :  = 


rjAD.eff  ’ 


(37) 
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Pressure-driven  convective  transport  is  described  sepa¬ 
rately  by  adding  a  term  to  the  diffusion  flux  densities  to  get 
the  overall  molar  flux  densities  tij  (mol  m-2  s-1),  analogue 
to  Eq.  (3): 


AD 


=  +  cfD<D 


(38) 


with  j  =  H20,  CH3OH,  C02. 

The  convective  velocity  nAD  (ms-1)  is  calculated  as  pre¬ 
sented  in  Section  3.2  (Eqs.  (4)  and  (5)).  As  the  carbon  dioxide 
and  methanol  concentrations  in  the  liquid  mixture  in  (AD)  are 
small  compared  to  the  water  concentration,  for  the  dynamic 
viscosity  as  function  of  local  temperature  and  pressure  a  cor¬ 
relation  for  pure  water  is  used  (calculation  see  Appendix  A. 5). 

The  binary  diffusion  coefficients  related  to  species-species 
interaction  were  determined  from  literature  correlations  (see 
Appendix  A. 9).  The  binary  diffusion  coefficients  related  to 
species-matrix  interaction  are  unknown.  But  for  liquid-phase 
transport  in  large  pores  (in  the  diffusion  layers  the  pore 
diameters  are  in  the  order  of  10-100  pm)  species-matrix 
interactions  are  small  compared  to  the  species-species 
interactions  for  diffusive  transport.  But  simply  skipping  the 
terms  with  the  species-matrix  binary  diffusion  coefficients 
in  Eq.  (36)  leads  to  numerical  problems,  as  the  resulting 
transport  matrix  [  Z?AD]  can  not  be  inverted  by  MatLab  due  to 
ill  conditioning  (too  close  to  singular).  Therefore,  the  values 
for  the  binary  species-matrix  diffusion  coefficients  were 
set  to  values  which  are  three  orders  of  magnitude  higher 
than  those  of  the  species-species  interactions.  Thus,  the 
numerical  problem  were  solved  while  the  influence  of  the 
wall  friction  on  the  individual  flux  densities  jf  becomes 
negligible. 

The  total  mass  flux  densities,  enthalpy  flux  densities  and 
conductive  heat  flux  densities  are  calculated  as  described  in 
Sections  3.2  and  3.3 


3.9.  Cathode  diffusion  layer  ( CD) 

Material  and  energy  balances  are  formulated  analogue  to 
(AD).  The  mobile  species  are  oxygen  (02),  nitrogen  (N2), 
water  (H20)  and  carbon  dioxide  (C02). 

The  major  difference  to  (AD)  results  from  the  fact,  that 
here  the  mixture  is  an  ideal  gas  and  not  a  liquid.  Therefore, 
no  total  mass  balance  is  needed  and  the  local  total  pressure 
pCD  js  gjven  by  the  sum  of  the  local  partial  pressures  p  f 

(Pa): 

p™  =  YjpT  =  RT™Y,cT  (39) 

j  j 


with  j  =  N2,  02,  H?0,  C02. 

For  diffusive  one-dimensional  mass  transport  in  a  mixture 
of  ideal  gases,  the  Maxwell-Stefan  Eq.  (2)  simplifies  to: 


1  dPj° 

RTcd  dz 


CD  -CD 


y,-  j 


yfDjfD 
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j-jCD.eff 
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CD  ,eff 


(40) 


jM 


with  j  =  N2,  02,  H20,  C02;  where  jJD  are  the  diffusive 
mass  flux  densities  in  (mol  m— 2s-1)  and  yf  are  the  gas 
mole  fractions. 

The  binary  diffusion  coefficients  can  be  easily  and  quite 
reliably  derived  from  several  correlations  (see  Appendix  A. 9) 
or  alternatively  from  the  kinetic  gas  theory.  The  species- 
matrix  diffusion  coefficients  are  calculated  using  the  equation 
for  Knudsen  diffusion: 


r-.CD,eff  _  Spores  4o?e  l%RTCD 

■>M  rCD  3  V 


(41) 


with  j  =  N2,  O2,  H2O,  CO2. 

In  Eq.  (41)  Mj  are  the  molecular  weights  (kg  mol-1)  and 
c/pjlr?.  (m)  is  the  mean  pore  diameter  in  the  matrix. 

The  mass  flux  densities  can  be  obtained  from  the  transport 
equations  following  the  same  numerical  method  as  described 
for  (AD).  Also  the  calculation  of  the  total  mass,  enthalpy  and 
heat  flux  densities  is  analogous.  Only  one  major  difference  to 
the  anode  side  has  to  be  accounted  for.  On  the  anode  side,  a 
liquid  mixture  (mainly  water)  is  pumped  through  the  channels 
of  (A).  Due  to  the  high  heat  capacity  of  water  and  the  rela¬ 
tively  high  heat  transfer  coefficients  between  the  liquid  mix¬ 
ture  and  the  channel  walls,  a  uniform  temperature  TA  can  be 
assumed.  On  the  cathode  side  the  situation  is  totally  different 
in  this  respect,  as  the  gas  mixture  has  a  heat  capacity  of  only 
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(42) 


which  is  negligible  compared  to  that  of  the  bipolar  plate 


(Pcp)  graphite 
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(43) 


The  effect  of  this  situation,  observed  in  the  experiments, 
is  that  the  air  pumped  through  the  cathode  bipolar  plate 
changes  its  temperature  only  marginally  (around  1-2  K  max¬ 
imum,  see  also  basic  model  assumptions)  and  remains  at  30- 
35  °C,  while  the  temperature  of  the  bipolar  plate  is  nearly 
the  same  as  that  on  the  anode  side  (up  to  90  °C  in  the  pre¬ 
sented  experiments).  Thermal  energy  is  transported  with  low 
resistance  through  the  planar  contact  areas  at  the  outer  gas¬ 
kets  and  through  the  solid  materials  of  the  MEA.  Therefore, 
on  the  cathode  side,  two  temperatures  have  to  be  accounted 
for  as  boundary/operating  conditions:  The  gas  temperature, 
which  nearly  equals  the  feed  temperature  Tc  (simplifyingly 
assumed  to  be  equal),  and  the  solid  temperature  7((alhon  •  For 
the  calculation  of  the  enthalpy  flux  (convective  heat  flux) 
between  (CD)  and  (C),  the  gas  temperatures  have  to  be  ac¬ 
counted  for,  as  the  gases  are  the  mobile  species.  For  the  cal¬ 
culation  of  the  conductive  heat  flux,  the  temperatures  of  the 
solid  matrices  have  to  be  used,  as  thermal  conduction  through 
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the  gas-phase  can  be  neglected  in  comparison  to  that  through 
the  solid  materials. 

But  this  distinction  is  only  necessary  for  the  energy  fluxes 
between  (CD)  and  (C).  Within  the  diffusion  layer  (CD)  still 
a  uniform  temperature  valid  for  gas  and  solid  matrix  can  be 
assumed,  as  here  the  flow  velocity  is  some  orders  of  magni¬ 
tude  lower  as  in  the  gas  channels  of  (C)  and  due  to  the  small 
pore  diameters  and  high  pore  tortuosity  an  intense  heat  ex¬ 
change  between  gas  and  solid  can  be  expected.  Therefore,  the 
temperature  can  also  be  expected  to  show  only  a  slight  slope 
within  (CD),  whereas  at  the  interface  between  (CD)  and  the 
gas  channels  in  (C)  it  can  change  dramatically,  depending  on 
the  cathode  and  anode  feed  temperatures. 


In  Eq.  (44)  for  the  polymer  material  the  lower  index  “M”  is 
used.  For  each  pair  of  species,  a  non-ideality  parameter  Xj.i  is 
required.  Crosslinking  of  the  polymer  material  is  accounted 
for  in  the  last  term  on  the  right  hand  side:  (Vm.cu  is  the  num¬ 
ber  of  sequential  single  polymer  chain  units  (i.e.  monomer 
units)  within  the  main  polymer  chain  between  two  cross¬ 
links,  Vm.cu  is  the  molar  volume  of  such  a  single  chain  unit 
in  (m3  mor1).  As  the  molar  volume  of  the  polymer  is  some 
orders  of  magnitude  higher  than  those  of  the  mobile  species, 
the  term  (1  —  Vj/Vf)  is  approximately  1.  In  the  following, 
three  species  are  accounted  for:  The  polymer  backbone,  wa¬ 
ter  and  methanol.  Therefore,  three  non-ideality  parameters 
are  needed,  which  were  determined  from  swelling  experi¬ 
ments  [15]  as 


3.10.  Membrane  electrode  assembly 


The  MEA  consists  of  the  polymer  electrolyte  membrane 
(M)  and  the  anode  and  cathode  catalyst  layers  (AD)  and 
(CD),  respectively.  These  three  layers  are  highly  inter¬ 
connected  by  the  proton-conducting  membrane  material 
content  within  the  catalyst  layers,  therefore  they  can  not  be 
treated  separately.  Strains  of  polymer  material  are  running 
through  the  catalyst  layers,  connecting  catalyst  particles  to 
the  membrane  on  the  ionic  conductor  level.  These  strains 
not  only  form  an  ionic  connection  for  mobile  protons  to  the 
membrane  (M),  but  also  with  respect  to  all  species  which  can 
enter  the  pores  within  this  material.  In  the  case  of  NAFION™ 
these  species  are  water  and  methanol,  whereas  the  membrane 
is  assumed  to  be  impermeable  for  all  gases  (carbon  dioxide, 
oxygen,  nitrogen).  This  assumption  is  justified  by  the  fact 
that  the  maximum  content  of  these  gases  in  the  pore  fluid 
of  the  membrane  is  limited  by  their  solubility  in  water  (and 
methanol),  which  is  very  low,  especially  at  high  operating 
temperatures. 

As  the  pore  volume  of  the  polymer-phase  within  the  cat¬ 
alyst  layers  is  small  compared  to  its  surface  open  to  the  free 
pores  within  the  catalyst  layers,  it  will  be  furthermore  as¬ 
sumed  that,  at  this  interface  between  the  catalyst  layers  and 
the  membrane,  phase  equilibrium  is  always  established. 

On  the  anode  side,  the  phase  equilibrium  can  be  described 
applying  a  UNIFAC  activity  model  [19]  for  the  liquid  in 
the  free  pores  of  (AC),  and  a  Flory-Huggins  activity  model 
[  1 7,20]  for  the  liquid  inside  the  polymer-phase  ( for  a  detailed 
description  see  [15]).  In  the  latter,  the  activity  of  a  species  j 
is  given  as  a  function  of  the  volume  fractions  s j  of  all  mo¬ 
bile  species  and  the  polymer  backbone,  treating  the  polymer 
backbone  and  the  mobile  species  as  a  liquid  mixture  (typical 
polymers  are  undercooled  liquids): 


aj  =  Sj  exp 


1 


Yi 


si  +  Xj.iSJ 


+ 


Vj  _ 

2  •  Am  ,CU  '  Em  ,CU 


(44) 


Xh2o,m  =  0.7177,  Xch3oh,m  =  0.1348  and  xh3o,ch3oh  =  1-3. 


The  protons,  strictly  speaking,  are  also  a  mobile  species, 
but  they  are  treated  separately,  as  will  be  presented  later. 

The  anode  phase  equilibrium  for  water  and  methanol, 
obtained  from  a  dynamic  mass  balance  model  using  both, 
UNIFAC  and  Flory-Huggins  activity  models,  can  be  ap¬ 
proximated  by  the  following  fitting  functions  for  (AC)  pore 
methanol  mole  fractions  below  0.03  (i.e.  typical  operation 
range  of  aDMFC)  [15]: 

£ch3oh  —  25.4831  ■  (-*ch3oh)3  +  4.2821  •  (xch3oh)“ 

+  1 .6354  •  Xch3oh  (45) 


4o  =  -104.9956  ■  (^3oh)3  +  20.9052  •  (*££30H)2 


■  2.6349  •  -Vch3oh 


0.4601 


(46) 


The  temperature  influence  was  found  to  be  negligible:  Ac¬ 
tivity  coefficients  predicted  from  the  UNIFAC  model  for  this 
system  show  only  a  weak  temperature  dependence,  and  in  the 
Flory-Huggins  model  temperature  influences  only  the  molar 
volumes,  which  are  also  nearly  independent  of  temperature 
for  liquids. 

A  second  phase  equilibrium  model  has  to  be  formulated 
for  the  interface  between  the  membrane  (M)  and  the  cath¬ 
ode  catalyst  layer  (CC)  pores.  Here  it  is  assumed,  that  in  the 
pores  of  (CC)  a  gas-phase  is  predominant  (i.e.  condensation 
of  water  is  neglected).  Experimental  data  are  available  from 
the  literature  (e.g.  [21,22])  for  the  equilibrium  relative  water 
content  of  NAFION™  (with  respect  to  the  number  of  fixed 
sulfonic  acid  groups  -SO3-): 


Am 


nm 

lyH2Q 

nr-so3~ 


(47) 


as  function  of  the  “water  vapour  activity”  (ideal  gas,  defined 
as  in  the  references  [21,22]): 


*  _  TH;U 

flH2o(g)  -  p^o(Ty 


(48) 
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It  is  possible  to  fit  a  polynomial  to  these  data  describing 
the  phase  equilibrium  here  [15]: 


^CCP  =  28.5  •  (off  -  0.35)3 


4H20(g) 

*,CC 

‘H20(g) 


+  5 -(4’^ -0-35) +  3 


(49) 


with 


*,cc  _  P  h2o 

H2°(g)  Pto(TCcy 


(50) 


Finally  some  attention  is  paid  to  viscous  flow  in  the  mem¬ 
brane  material.  As  the  pores  within  NAFION™  are  extremely 
small  (in  the  range  of  a  few  nanometers)  the  question  arises 
whether  small  pressure  differences  over  the  DMFC  (typically 
maximum  2  bars)  can  lead  to  a  significant  viscous  flow  con¬ 
tribution.  This  question  can  be  addressed  by  some  simple 
calculations  [15].  The  outcome  is,  that  a  pressure  difference 
of  roughly  A pM  ~  25  bars  would  be  required  to  obtain  water 
crossover  fluxes  of  the  order  of  magnitude,  which  was  ob¬ 
served  in  the  experiments.  As  in  real  DMFC  operation  the 
maximum  pressure  differences  do  not  exceed  1-2  bars,  one 
can  conclude  that  viscous  flow  can  only  contribute  a  few  per¬ 
cent  to  the  measured  membrane  crossover  fluxes.  Therefore, 
it  is  justified  to  neglect  pressure-driven  flow  in  the  membrane 
model. 


3.10.1.  Anode  catalyst  layer  (AC) 

In  the  anode  catalyst  layer  one  finds  several  phases  which 
are  highly  interconnected:  free  pores,  polymer-phase  and 
electron  conductor.  Water,  methanol  and  carbon  dioxide  are 
the  mobile  species  within  the  free  pores.  The  catalyst  layers 
are  modelled  as  concentrated  parameter  systems.  The  species 
mass  balances  of  (AC)  are  given  by: 


dcf 

df 


nju(zAU  =  dAU )  -  nf(zM  =  0) 

oAC  ,/AC 
c  pores 


+  VaJ  ra 


(51) 


with  j  —  FbO,  CH3OH,  CO2  and  the  stoichiometric  coeffi¬ 

cients  of  the  anodic  electrochemical  methanol  oxidation: 


fa.CHjOH  =  —  1  Va,H20  =  —  1 

va,C02  =  +1 * * * *  va  ,H+  =  +6  (52) 


The  rate  of  the  electrochemical  methanol  oxidation  is  for¬ 
mulated  as  a  Butler- Volmer  equation: 


ra  =  ka 


.’All, oil  exP 


/  aa6F  \ 
V  RTAC 


-  xc52  exp 


/  (1  —  «a)6  F  \ 
V  RTAC  11'd) 


(53) 


In  Eq.  (53)  kd  is  the  rate  constant  (molm-3  s_1)  (value 
for  simulation:  kd  =  6  x  10-3  mol  m-3  s_1),  aa  is  the  charge 
transfer  coefficient  (set  to  aa  =  0.1)  and  r/a  is  the  anode  over¬ 
potential  (V).  The  water  mole  fraction  is  not  included  in  Eq. 
(53)  due  to  the  fact  that  it  can  be  assumed  to  be  close  to  unity 


and  is  not  changing  significantly.  Due  to  the  very  complex 
reaction  mechanism  this  rate  equation  can  of  course  be  only  a 
first  approach.  More  realistic  models  of  the  reaction  kinetics 
have  to  account  for  methanol  oxidation  reaction  intermedi¬ 
ates  and  adsorption  and  desorption  phenomena  on  the  binary 
anode  catalyst.  Such  more  detailed  kinetic  models  were  e.g. 
proposed  in  [  1 ,2] ,  but  the  number  of  free  model  parameters  to 
be  determined  is  considerable,  as  is  the  increase  in  required 
computation  time. 

The  concentrations  of  water  and  methanol  in  the  polymer- 
phase  within  (A)  (denoted  as  ACP)  are  calculated  assuming 
phase  equilibrium  with  the  pores  in  (AC),  as  described  above 
(Section  3.10). 

The  total  mass  balance  in  (AC)  is  given  by 

d  PAC  ...  ^sound  ( M c\\  ^„AD 

df  =  “  dKC  =  0)  —  mtot(z  =d  )j 

(54) 

As  the  cell  is  operated  galvanostatically  (i.e.  the  cell  cur¬ 
rent  density  fceii  is  a  known  operating  parameter),  the  charge 
balance  can  be  formulated  quasi-stationary: 

0  =  /cell  -  iM(zM  =  0).  (55) 

Due  to  the  quasi-stationarity  of  the  charge  balance,  the 
mass  and  charge  balances  are  not  coupled.  Therefore,  the 
current  density  of  the  anodic  reaction,  ia,  expressed  in  terms 
of  the  anodic  rate  expression,  Eq.  (53),  is  identical  to  the 
known  cell  current  density: 

i-d  —  d  ■  £p0res  -6  ■  F  ■  r a(/^a)  =  /cell-  (56) 

From  Eq.  (56)  and  the  anodic  rate  expression,  Eq.  (53), 
the  anode  overpotential  can  be  calculated  numerically  by  re¬ 
cursion. 

Finally,  the  energy  flux  densities  at  the  interfaces  between 
(AC)  and  (AD)  and  (AC)  and  (M),  respectively,  are  accounted 
for  in  the  energy  balance  of  (AC): 

(( eM(zM  =  0)  -  eAD(zAD  =  dAD)) 
d tAC  _  +  (qM(zM  =  0)  -  qAD(zAD  =  dAD 

dr  (pFp)ACdAC 

3.10.2.  Cathode  catalyst  layer  ( CC) 

Very  similar  to  (AC)  also  in  (CC)  spatially  concentrated 
balances  for  mass,  charge  and  energy  are  formulated.  The 
mass  balances  for  oxygen,  nitrogen,  water  vapour  and  carbon 
dioxide  in  terms  of  molar  concentrations  are: 

dcf  „M(Z M  =  dM )  -  nf(zM  =  0) 

fovcsdCC 

“I"  Vc,jfc  “I"  Across, j  Across  (58) 

with  j  =  N2,  O2,  H2O,  CO2 

In  Eq.  (58)  not  only  the  desired  electrochemical  oxygen 
reduction  has  to  be  accounted  for,  but  also  the  undesired  direct 
methanol  oxidation: 
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(a)  Oxygen  reduction  (consumes  H+  and  e  ): 

1.5  02  +  6H+  +  6e“  -*  3H20  (59) 

(b)  Methanol  oxidation  (produces  H+  and  e_): 

CH3OH  +  0.5O2  -»  C02  +  4H+  +  4e“  (60) 

The  stoichiometric  coefficients  of  these  reactions  are 


As  at  the  cathode  two  electrochemical  reactions  take  place 
simultaneously  (oxygen  reduction  with  the  charge  production 
rate  rc,  and  methanol  oxidation  with  the  charge  production 
rate  rcross),  a  charge  balance  has  to  be  formulated  either  for 
the  electron  or  the  proton  conductor-phase.  For  the  electron 
conductor  phase  this  quasi-stationary  charge  balance  is  given 
by: 


fc,o2  =  — 1-5  Uc,n2  =  0  vc,h2o  =  +3 

vc,co2  =  0  Vc  H+  =  -6 

and 


Across,  02  —  0.5  Vcross,N2  —  0  VCross,H20 

Across,  C02  =  +1  l’cross,H+  =  +4- 


0 


It  is  known  that  the  reaction  mechanism  of  the  electro¬ 
chemical  oxygen  reduction  is  complicated  with  several  pos¬ 
sible  intermediates.  However,  as  in  (AC),  in  order  to  keep  the 
number  of  unknown  parameters  small,  again  a  Butler- Volmer 
type  equation  is  applied: 


rc  =  kc 


Pen  " 

105  Pa 


1.5 


exp 


f  ac  6F 
\Rfcc 


(  (l-ac)6  F  \ 
-  6XP  RTCC  * ) 


(61) 


0  —  ('cell  +  (('cross  T  ("c)-  (65) 

In  Eq.  (65)  the  current  density  of  the  cathode  oxygen  re¬ 
duction  reaction  is  given  by 

ic  =  dccs%s6Frc  (66) 

(6  exchanged  electrons  per  net  reaction)  and  the  current  den¬ 
sity  of  the  oxidation  of  crossover  methanol  analogous  by 

Cross  =  d  £pores^^"cross-  (67) 

(4  exchanged  electrons  per  net  reaction).  From  Eqs.  (65)- 
(67)  together  with  the  cathode  rate  expression,  Eq.  (61),  the 
cathode  overpotential  i]c  can  be  determined  numerically  in  a 
recursion,  similar  to  the  anode  overpotential. 

Finally,  analogue  to  (AC),  the  energy  balance  for  (CC)  is 
formulated  as: 

((eCD(„CD  =  0)  _  gM(„M  =  dMy 

d tcc  +  (qCD(zCD  =  0)  -  t/M(2M  =  r/M))) 

~dT  ~  (pcp)ccdcc  ‘  (68) 


In  Eq.  (61)  kc  is  the  rate  constant  (molm-3  s-1)  (value 
for  simulation:  kc  =  1.27  x  10~21  molm-3  s-1),  ac  is  the 
charge  transfer  coefficient  (set  to  ac  =0.5)  and  rjc  is  the  cath¬ 
ode  overpotential  (V). 

Methanol  is  assumed  to  be  immediately  consumed  when 
coming  into  contact  with  the  cathode.  Therefore,  its  concen¬ 
tration  in  (CCP)  and  (CC)  drops  to  zero.  Under  these  con¬ 
ditions  the  rate  of  the  direct  oxidation  of  methanol  at  the 
cathode,  rcross  is  proportional  to  the  methanol  flux  from  (M) 
to  (CC): 


r cross 


n 


M 

ch3oh 


(zM  =  dm) 


IVK 


rlCCpCC 
“  c  pores 


(62) 


The  concentration  of  water  in  (CCP)  is  calculated  from 
the  equilibrium  condition  presented  in  Section  3.10. 

The  overall  pressure  pcc  is  calculated  from  the  concen¬ 
trations  of  all  four  gas  species  <f  according  to  the  ideal  gas 
law: 


pcc  =  RTccJ2fC 

j 


(63) 


Following  the  same  argumentation  as  for  (AC),  the  total 
charge  balance  is  formulated  quasi-stationary,  and  is  decou¬ 
pled  from  the  mass  balances  as  the  cell  is  operated  galvano- 
statically: 

0  =  (Ceii  -  CM(ZM  =  dM )  (64) 


3.10.3.  Polymer  electrolyte  membrane  (M) 

The  PEM  (M)  is  a  one-dimensional  transport  element  like 
the  diffusion  layers  (AD)  and  (CD).  The  mass  balances  for 
water  and  methanol  are  given  by: 


dcf  dnf 

3 1  dz 


(69) 


with  j  =  H20,  CH3OH. 

But  the  material  structure  and  the  occurring  physical  phe¬ 
nomena  are  much  more  complex.  The  PEM  has  not  a  con¬ 
stant  porosity,  but  one  that  strongly  depends  on  the  local 
water  and  methanol  content.  The  relative  water  content  A, 
Eq.  (47),  can  have  values  between  A  —  0  (totally  dry  mem¬ 
brane)  and  A  &  30  (fully  swollen  with  water  and  methanol  at 
room  temperature),  depending  on  temperature  and  other  con¬ 
ditions.  Different  water  and  methanol  contents  result  in  dif¬ 
ferent  degrees  of  swelling,  and  therefore  different  thicknesses 
and  porosities.  Therefore,  it  is  not  suitable  to  formulate  mass 
balances  in  molar  concentrations,  as  these  refer  to  a  constant 
overall  volume.  It  is  more  convenient  to  use  a  concentration 
measure  which  refers  to  the  constant  cross-sectional  area  of 
the  cell,  As  (m2).  This  molar  density  N j  (mol  m-2)  is  defined 
as  quotient  of  the  total  molar  amount  of  species  j,  Nj  (mol), 
and  the  cell  cross-sectional  area: 


(70) 
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For  the  mass  transport  equations  (see  below)  a  variety 
of  other  concentration  measures  are  needed.  The  necessary 
conversions  are  given  in  Appendix  E. 

Using  this  introduced  molar  density,  the  mass  balances  for 
the  control  volumes  (index  k)  in  the  discretised  model  have 
the  form 


rl  P 

_  M 
- 1 -  —  n  i  k 

d  t  h 


'  nj,k+l 


(71) 


with  j  =  H20,  CH3OH. 

In  Eq.  (71),  swelling  of  the  membrane  is  assumed  to  be  in 
the  steady  state.  A  total  mass  balance  is  not  formulated,  as 
pressures  are  not  discussed  within  the  membrane. 

The  charge  balance  is  considered  in  steady  state  (quasi¬ 
stationary): 


0=  - 


djM 

3  z 


(72) 


As  the  mass  balances,  the  energy  balance  is  quite  similar 
to  that  in  the  diffusion  layers.  The  only  difference  is  that  it  has 
to  account  not  only  for  convective  and  conductive  heat  fluxes, 
but  also  for  Joule  heating  resulting  from  the  electric  resistance 
of  the  membrane  caused  by  friction  between  mobile  charged 
species  and  the  immobilised  counter  charges  at  the  pore  walls. 
The  heat  production  due  to  Joule  heating,  £^ule  (J  m-3  s_l),  is 
proportional  to  the  electric  current  density  iM  and  the  gradient 
of  the  electric  potential  </>M: 


^Joule 


•M 


3  4>m 

3  z 


With  all  this,  the  energy  balance  finally  is 

3 TM  _  1  T  deM  3 qM  ,m3 <pM' 

3 1  (pcp)M  dz  3 z  dz 


(73) 


(74) 


In  the  following,  transport  equations  for  mass  and  energy 
are  formulated  for  (M).  The  charge  transport  is  expressed  in 
terms  of  the  proton  flux.  Mass  transport  is  described  using 
the  Maxwell-Stefan  approach,  as  has  already  been  done 
for  the  diffusion  layers.  But  in  the  PEM,  a  more  complex 
formulation  of  the  driving  forces  has  to  be  chosen.  First, 
the  migration  term  has  to  be  included,  as  one  of  the  mobile 
species  (protons)  is  charged  and  an  electric  field  is  present 
within  the  membrane.  Second,  the  diffusive  term  has  to 
account  for  the  highly  non-ideal  behaviour  of  the  mobile 


species  within  the  membrane  pores,  i.e.  the  gradients  of 
the  chemical  potentials  have  to  be  used  as  driving  force, 
which  are  equal  to  the  gradients  in  the  species  activities 
(calculated  using  the  presented  Flory-Huggins  activity 
model).  The  pressure-dependency  of  the  chemical  potentials 
as  well  as  viscous  flow  due  to  pressure  differences  across  the 
membrane  can  be  neglected  as  was  shown  in  Section  3.10. 

All  this  leaves  the  following  form  of  the  Maxwell-Stefan 
equations  for  the  mobile  species  (j  —  H+,  H2O,  CH3OH): 
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In  Eq.  (75),  as  superficial  viscous  flow  due  to  pressure 
gradients  is  not  accounted  for,  the  total  molar  flux  densities 
Hi  appear  in  the  friction  terms  on  the  right  hand  side  of  the 
equation.  Three  of  the  six  binary  diffusion  coefficients  D^'eff 
are  taken  from  the  literature  [23]  as  a  starting  point,  using  the 
general  expression 


Dff(T,  Am)  =  D f;-AM  exp 


fl 

j_y 

R 

Kt 

Te  J 

(76) 


with  the  reference  temperature  Te  —  298  K.  The  values  of 
the  standard  diffusion  coefficients  and  the  activation  energies 
are  given  in  Table  2.  The  two  additional  binary  diffusion  co¬ 
efficients  for  methanol/pore  wall  and  methanol/protons  are 
formulated  in  the  same  way,  using  the  literature  values  for 
water/pore  wall  and  water/protons  as  a  first  estimate,  respec¬ 
tively.  The  coefficient  for  methanol/water  is  calculated  using 
the  free  solution  correlation  of  Hayduk  and  Minhas  for  so¬ 
lutes  in  aqueous  solutions  [34d] : 


E>ch3oh,h2o  ~  h2o 

=  1-25  -  10-12(Uch3oh  -  0.292)r152()7™or  (77) 

9.58 

with  e  =  - - 1.12. 

VcihOii 

All  details  on  derivation  of  the  binary  diffusion  coefficients 
are  given  in  Appendix  A. 9. 


Table  2 

Parameters  for  calculation  of  binary  diffusion  coefficients  in  Nafion™  (M  represents  solid  matrix/pore  wall) 


Species  pair  ( i / j) 

Original  (*  =[23])  parameters 

Adjusted  parameters 

Df  .  (m2  s  1 ) 

£5  (kJ  mol-1) 

Df.  (m2  s_1) 

£5  (kJ  mol  1 ) 

h2o/h+ 

0.85  x  10~10  (*) 

10.54  (*) 

0.15  x  10~10 

10.54 

h2o/m 

0.55  x  10~u  (*) 

20.25  (*) 

0.20  x  10~n 

50.25 

H+/M 

0.22  x  10-10  (*) 

10.54  (*) 

0.22  x  10~10 

10.54 

ch3oh/h+ 

Identical  to  ThO/tU 

0.60  x  10~10 

8.43 

ch3oh/m 

Identical  to  H2O/M 

5.00  x  lO^11 

25.13 

ch3oh/h2o 

1.25  x  10-12 

- 

5.00  x  10-12 

- 
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As  the  flux  density  of  protons  is  given  by  the  electric 
current  density  iceii  using  Faraday’s  law 

:M 
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,?H+  -  Z7  - 


*cell 


(78) 


only  the  flux  densities  of  water  and  methanol  have  to  be  de¬ 
termined.  In  this  case,  the  flux-implicit  transport  equations, 
Eq.  (75),  can  be  easily  transformed  into  a  flux-explicit  form 
by  rearranging: 
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The  electric  potential  gradient  d(f>M/dz  in  the  mem¬ 
brane  material,  due  to  the  transport  resistance  to  the  proton 
flux  (“Ohmic  drop”  over  membrane)  is  obtained  from  the 
Maxwell-Stefan  equation  for  the  protons  as: 
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The  problem  with  Eq.  (89)  is,  that  for  the  protons  no  Flory- 
Huggins  non-ideality  parameters  are  known.  Therefore,  in  the 
above  equations,  the  activity  of  protons  is  approximated  by 
the  mole  fraction  of  protons  in  the  pore  liquid,  x^+ . 

For  the  energy  balance  in  (M),  transport  equations  are 
needed  for  thermal  conduction  and  convective  heat  transport. 
Both  are  similar  to  the  equations  for  (AD)  and  (CD): 

d  tm 

Conductive  heat  flux  :  qM  =  -lM'eff - ,  (94) 

dz 

Enthalpy  flux  (convective  heat  flux)  :  eM  —  i  f  h  ;(TM). 

j 

(95) 


4.  Simulation  results 

The  presented  DMFC  model  was  implemented  in  MatLab 
using  the  solver  odel5s  to  carry  out  dynamic  simulations  into 
steady  state  at  a  variety  of  operating  conditions.  The  param¬ 
eters  which  were  varied  are  the  anode  feed  temperature  TAF 
and  the  cell  current  density  /ceii-  Fig.  4  presents  some  of  the 
results  of  these  steady  state  simulations  and  corresponding 
experimental  data  for  anode  feed  temperatures  from  30  up 
to  90  °C  at  varying  cell  current  densities.  All  other  operating 
conditions  are  given  in  the  figure.  It  has  to  be  emphasised  that 
all  simulation  results  are  obtained  using  the  same  set  of  pa¬ 
rameters  (presented  in  the  preceding  sections  and  the  respec¬ 
tive  appendices).  The  parameters  for  the  calculation  of  the 
binary  diffusion  coefficients  are  based  on  those  found  in  the 
literature  [23],  but  have  been  slightly  adjusted  in  order  to  get 
a  better  fit  of  the  experimental  data  (values  given  in  Table  2). 
The  activation  energies  used  (right  column  in  Table  2)  are  in 
the  typical  range  for  diffusive  transport.  Interesting  is  the  high 
value  for  the  pair  water/membrane.  Here  possibly  additional 
thermal  effects  are  reflected,  like  those  related  with  solvation. 

As  one  can  see  from  Fig.  4,  in  general  a  reasonable  ap¬ 
proximation  to  the  experimental  steady  state  results  has  been 
achieved.  The  simulation  results  for  the  membrane  crossover 
flux  densities  as  well  as  for  the  current-voltage  curves  are 
in  the  orders  of  the  experimental  data.  Also  the  trends  are 
predicted  correctly,  i.e.  water  crossover  fluxes  increase  with 
current  density  and  methanol  crossover  fluxes  decrease  with 
current  density.  Especially  for  the  methanol  crossover,  all 
simulation  results  are  within  or  close  to  the  error  bars  of  the 
experimental  data  (error  bars  based  on  evaluation  of  the  ac¬ 
curacy  of  the  sensors  used  to  determine  the  crossover  fluxes 
[15]).  As  the  methanol  crossover  plays  a  key  role  for  the 
performance  of  the  DMFC,  its  correct  prediction  is  one  of 
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Fig.  4.  Experimental  (symbols)  and  simulated  results  (lines)  for  a  single-cell  DMFC.  Left:  crossover  flux  densities  (left  y-axis:  H2O,  right  y-axis: 
CH3OH);  right:  current-voltage-curves  (rAF  =  30,  . . . ,  90  °C,  pA  =  pc  =  1.7  bar,  c^OH  =  1  mol  dm-3,  c^q2  =  1  |xmoldm-3,  Faf  =  0.5  dm3  min-1, 
Fcf  =  0.5  scbm  h-1 ,  cathode  feed:  air  with  dew  point  3  °C  at  1  bar). 


the  most  important  benchmarks  of  any  mathematical  DMFC 
process  model. 

Significant  deviations  from  the  experimental  observations 
exist  for  the  crossover  water  flux  densities  at  very  low  and 


very  high  temperatures.  The  experimental  crossover  water 
fluxes  through  the  membrane  show  an  increasing  gradient 
with  increasing  current  density.  Such  a  behaviour  can  not 
be  explained  by  the  Maxwell-Stefan  model,  as  the  model 
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allows  an  utmost  only  linear  dependency  between  flux  and 
current  density.  Also  the  binary  diffusion  coefficients  within 
the  membrane  material  are  formulated  such,  that  they  are 
increasing  linearly  with  the  local  relative  water  content  (a 
very  recent  publication  supports  this  assumption  [24]).  Fur¬ 
ther  optimisation  of  the  parameters  for  the  calculation  of  the 
binary  diffusion  coefficients  can  be  expected  to  yield  even 
better  approximations  to  the  experimental  crossover  fluxes. 
For  this  purpose,  numerical  optimisation  methods  can  be  ap¬ 
plied.  But  for  this  task,  first  the  computational  time  of  the 
model  should  be  reduced  by  either  optimising  the  source 
code,  or  implementing  the  model  in  other,  faster  solver  tools 
than  MatLab.  Finally,  also  the  applied  Flory-Huggins  activ¬ 
ity  model  influences  the  simulation  results  in  terms  of  the 
membrane  crossover  fluxes.  The  Flory-Huggins  model  was 
originally  formulated  for  mixtures  of  uncharged  polymers 
and  uncharged  solvents,  therefore  its  application  to  a  poly¬ 
mer  electrolyte  demands  further  discussion. 

After  analysing  the  simulated  crossover  fluxes,  now 
the  simulated  current-voltage  curves,  overpotentials  and 
potentials  drops  over  the  membrane  are  discussed.  The 
results  are  presented  in  the  left  column  of  Fig.  4.  Generally 
it  can  be  observed  that  the  cathode  overpotentials  are  very 
high,  even  for  open  circuit  conditions.  This  is  due  to  the 
oxidation  of  crossover  methanol  at  the  cathode,  and  the 
resulting  mixed  potential  formation.  The  influence  of  the  cell 
current  density  on  the  cathode  overpotential  is  only  small.  In 
contrast  to  this,  the  anode  overpotential  varies  significantly 
with  the  cell  current  density.  At  open  circuit  conditions, 
the  anode  overpotential  is  close  to  zero,  for  high  current 
densities,  values  around  0.25  V  are  found.  Finally,  the  total 
potential  drop  over  the  membrane  shows  a  nearly  linear 
increase  with  only  a  small  slope.  What  seems  remarkable  is 
the  fact,  that  it  is  not  zero  at  open  circuit  conditions.  This 
can  be  explained  with  the  diffusive  flux  of  water  through  the 
membrane,  which  also  takes  place  when  there  is  no  electrical 
current  flow.  The  water  molecules  exert  a  drag  on  the  protons 
in  the  membrane,  but  the  protons  are  withheld  by  electro¬ 
static  forces  between  them  and  their  counter-ions  bound 
to  the  membrane  material.  An  electric  field  is  produced  by 
this  phenomenon,  which  is  often  referred  to  as  streaming 
potential. 

Another  observation  from  Fig.  4  is,  that  the  experimental 
open  circuit  cell  voltages  increase  with  the  cell  temperature, 
while  the  model  predicts  a  decrease.  This  is  to  a  large  extent 
due  to  the  fact,  that  in  the  model  the  open  circuit  cell  voltage 
is  calculated  from  thermodynamics  using  some  simplifying 
assumptions  [15].  These  thermodynamic  relations  exhibit  a 
decrease  of  the  open  circuit  cell  voltage  with  increasing  tem¬ 
perature.  The  difference  between  the  thermodynamic  and  the 
real  behaviour  results  from  the  fact,  that  in  reality  not  a  one- 
step  total  methanol  oxidation  takes  place  at  the  anode  (as  is 
assumed  in  the  thermodynamic  considerations),  but  a  very 
complex,  multi-step  network  of  adsorption  and  desorption 
processes  and  reaction  intermediates  exists.  A  better  predic¬ 
tion  of  the  open  circuit  voltage,  based  on  a  more  realistic 


situation  at  both  electrodes,  would  significantly  enhance  the 
prediction  of  the  current-voltage  curves. 

Nonetheless,  for  moderate  cell  current  densities,  the  model 
predicts  slopes  of  the  current-voltage  curves,  which  are  close 
to  the  experimental  results.  In  this  regime,  the  cell  behaviour 
is  dominated  by  mass  transfer  phenomena  within  the  mem¬ 
brane,  which  seems  to  be  reasonably  represented  by  the 
model. 

At  high  cell  current  densities,  finally,  the  predicted  cell 
voltages  are  much  higher  than  the  observed  experimental  val¬ 
ues.  Also,  the  experimental  results  show  a  typical  limiting 
current  behaviour  (breakdown  of  the  cell  voltage),  while  the 
model  shows  such  limiting  current  behaviour  only  for  signif¬ 
icantly  higher  cell  current  densities  (not  shown  in  the  plots  in 
Fig.  4).  Here  it  becomes  evident  that  the  model  is  based  on 
severe  simplifications  with  respect  to  mass  transport  in  both 
diffusion  layers.  The  model  does  not  account  for  the  possi¬ 
ble  coexistence  of  two  phases  (gas  and  liquid)  within  both 
diffusion  and  catalyst  layers,  although  it  is  well-known  from 
various  experimental  observations.  Carbon  dioxide  bubbles 
are  released  from  the  anode  diffusion  layer  [25,26],  conden¬ 
sation  of  water  can  occur  inside  the  cathode  pore  structure 
(so-called  cathode  flooding)  [8].  Both  phenomena  lead  to  in¬ 
creased  transport  resistances  for  the  fuel  (methanol)  and  ox¬ 
idant  (oxygen)  to  the  respective  electrodes  and  they  are  both 
most  important  for  high  current  densities,  i.e.  for  the  limit¬ 
ing  current  behaviour.  Such  two-phase  transport  behaviour 
has  therefore  to  be  included  in  a  DMFC  model  if  a  realistic 
simulation  of  the  limiting  current  behaviour  is  to  be  achieved. 

As  the  model  is  one-dimensional  perpendicular  to  the  cell 
plane,  profiles  through  the  DMFC  are  obtained  for  concen¬ 
trations,  temperature,  pressure  and  all  presented  fluxes.  Ex- 
emplarily  selected  steady  state  profiles  are  presented  in  Fig. 
5  for  an  anode  feed  temperature  of  60  °C  and  a  cell  current 
density  of  200mA cm-2.  All  other  parameters  are  the  same 
as  those  given  in  Fig.  4.  In  Fig.  5  the  ordinates  show  the  real 
cell  geometry  with  respect  to  the  thicknesses  of  the  different 
layers  of  the  DMFC.  The  vertical  lines  represent  the  limits  of 
the  control  volumes,  illustrating  the  spatial  discretisation  of 
both  diffusion  layers  (AD,  CD)  and  the  membrane  (M).  One 
can  see  that  both  diffusion  layers  are  represented  by  five  con¬ 
trol  volumes  each,  and  that  the  membrane  is  discretised  into 
ten  control  volumes.  It  is  also  apparent  that  the  thicknesses 
(i.e.  the  volumes)  of  the  diffusion  layer  control  volumes  are 
constant,  as  these  layers  consist  of  a  rigid  solid  matrix.  In 
contrast  to  this,  the  thicknesses  (and  therefore  also  the  total 
volumes)  of  the  membrane  control  volumes  change  due  to  dif¬ 
ferent  water  and  methanol  contents,  representing  the  swelling 
behaviour  of  the  membrane  material.  Simplifyingly,  in  the 
model  all  volume  changes  due  to  swelling  only  influence  the 
thicknesses  of  the  control  volumes  along  the  model  coordi¬ 
nate  perpendicular  to  the  cell  plane.  It  has  to  be  mentioned 
that  the  thicknesses  of  the  anode  and  cathode  channels,  (A) 
and  (C),  respectively,  do  not  represent  their  real  dimensions. 

One  can  see  that  for  the  diffusion  layers, (AD)  and  (CD), 
nearly  linear  concentration  and  partial  pressure  profiles  are 


pressure,  p /  [Pa]  part,  press,  pNs,  Pq2,  PH2<),  PCoAPa1  activities, aH+,  10.aCHsOH/[-]  pore  conc.C  chrOH-  Cc02  /[mol/m3] 
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Fig.  5.  Simulated  steady  state  profiles  (T’AF  =  60  °C,  /Ceii  =  200  mA cm  2,  all  other  operating  conditions  as  in  Fig.  4).  Ordinate  represents  real  cell  geometry, 
vertical  lines  are  limits  of  control  volumes  (abbreviations  A,  AD,  etc.  used  to  denote  the  DMFC  layers,  see  list  of  symbols). 
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obtained.  The  slopes  of  the  partial  pressure  profiles  in  (CD) 
are  only  small.  Also  the  total  pressure  differences  between 
the  supply  channels  and  the  respective  catalyst  layers  are 
only  in  the  order  of  a  few  mPa.  Obviously,  according  to  the 
here  applied  model,  mass  transport  in  the  diffusion  layers 
is  quite  fast,  especially  on  the  cathode  (i.e.  gas)  side  (as  was 
already  pointed  out  above  referring  to  the  limiting  current 
behaviour). 

Shown  in  Fig.  5  are  also  the  phase  equilibria  for  water  and 
methanol  within  the  catalyst  layers  between  the  free  pore 
concentrations  and  the  equilibrium  concentrations  within  the 
membrane  material  (as  described  in  chapter  3).  One  can  see 
that  on  the  anode  side  the  methanol  concentration  in  the  mem¬ 
brane  pores  is  slightly  higher  than  that  in  the  free  pores,  while 
the  water  concentration  in  the  membrane  pores  is  a  little  lower 
than  that  in  the  free  pores  (see  circles  in  the  upper  two  plots  of 
Fig.  5).  The  phase  equilibrium  for  water  in  the  cathode  cata¬ 
lyst  layer  (CC)  is  also  highlighted  by  circles  in  the  respective 
plots. 

The  most  interesting  concentration  profiles  develop  within 
the  membrane  (dashed  concentration  profiles  in  the  upper  two 
diagrams  of  Fig.  5).  The  methanol  pore  concentration  shows 
a  strongly  bent  profile  in  the  direction  of  the  overall  flow,  i.e. 
towards  the  cathode.  This  makes  sense  as  methanol  is  dragged 
along  with  the  water  flow  (diffusion  and  electro-osmosis). 
Also  the  water  profile  is  slightly  bent  in  the  same  manner  due 
to  electro-osmotic  transport.  Nonetheless,  diffusion  remains 
the  major  mode  of  transport  for  methanol  and  water.  Most 
interesting  is  the  big  difference  in  the  water  content  between 
anode  and  cathode  side  of  the  membrane.  While  on  the 
anode  side  a  relative  water  content  of  around  26  is  reached, 
on  the  cathode  side  only  values  around  4  are  found.  This  is 
due  to  the  operation  of  the  cell  with  dry  air  at  high  flow  rates. 
Water  is  transported  away  from  the  cathode  catalyst  layer 
(CC)  very  efficiently,  drying  out  this  side  of  the  membrane 
according  to  the  phase  equilibrium  relation  (Eq.  (49)).  This 
change  in  water  content  is  also  illustrated  by  the  decreasing 


thickness  of  the  membrane  control  volumes  from  anode 
to  cathode. 

The  conductivity  of  the  membrane  is  given  by  the  friction 
exerted  on  the  moving  protons.  This  friction  is  represented 
by  the  binary  diffusion  coefficients,  which  in  turn  are  func¬ 
tions  of  the  local  water  content.  Therefore,  also  the  proton 
conductivity  is  a  function  of  the  local  water  content  and  thus 
varies  locally.  The  same  is  true  for  the  potential  gradient  in 
the  membrane. 

Finally,  the  temperature  profile  exhibits  only  very  small 
gradients  over  the  inner  layers  of  the  DMFC.  The  total  tem¬ 
perature  difference  between  anode  channels  and  outer  side  of 
the  cathode  diffusion  layer  (CD)  is  less  than  3  K.  Only  the  air 
in  the  cathode  channels  has  a  much  lower  temperature  close 
to  its  inlet  temperature,  due  to  the  short  residence  time  and 
the  small  heat  exchange  coefficients  between  channel  walls 
and  gas  (see  discussion  in  Section  3.9). 

From  the  simulation  results,  it  is  possible  to  evaluate 
the  importance  of  the  different  mass  transport  contributions 
(driving  forces  and  friction)  in  the  generalised  Maxwell- 
Stefan  framework,  Eq.  (1).  Table  3  presents  the  quintessence 
of  this  evaluation.  In  the  top  line,  the  complete  Eq.  (1)  is 
given.  In  the  following  rows  the  importance  of  the  individual 
terms  of  the  generalised  Maxwell-Stefan  equation  is  indi¬ 
cated  for  each  of  the  three  mass  transport  related  layers  of  the 
DMFC  and  each  mobile  species  by  “++"  (very  important), 
“+”  (moderately  important)  and  blanks  (not  important/ 
negligible). 

Obviously,  multi-component  diffusion  represented  by  the 
gradient  in  the  chemical  potentials  as  driving  force,  left  term, 
and  both  friction  terms  (species-species  and  species-matrix), 
right  side  of  Eq.  (1),  are  the  most  important  influencing  factors 
for  mass  transport  of  all  mobile  species. 

The  pressure-dependence  of  the  chemical  potentials 
(second  driving  force  term  on  the  left  hand  side  of  Eq.  (1)) 
is  negligible  in  all  DMFC  layers.  This  is  generally  justified 
for  liquid  phases  if  no  large  pressure  gradients  exist.  This 


Table  3 


Importance  of  mass  transport  contributions  in  Eq.  (1)  (driving  forces  and  friction  terms):  ++  =  very  important,  +  =  important,  blanks=  not  important/negligible 
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term  is  only  relevant  for  applications  with  extremely  high 
pressure  differences  as  they  can  be  found  e.g.  in  reverse 
osmosis  and  pervaporation  processes.  On  the  gas  side  of 
the  presented  model,  i.e.  in  the  cathode  diffusion  layer 
(CD),  it  is  also  negligible  due  to  the  very  low  pressure 
gradient. 

The  third  driving  force  is  pressure-driven  convection.  This 
is  accounted  for  in  both  diffusion  layers  and  has  shown  to  play 
an  important  role  especially  for  the  excess  components  in  the 
respective  fluid  mixtures,  i.e.  water  in  (AD)  and  nitrogen  in 
(CD).  Within  the  polymer  electrolyte  membrane  this  term  is 
neglected  due  to  the  low  hydraulic  permeability. 

Finally,  the  electric  field  as  driving  force  only  applies  to 
protons  as  the  only  mobile  charged  species. 


5.  Conclusions 

Based  on  a  systematic  approach,  a  one-dimensional  pro¬ 
cess  model  of  a  DMFC  has  been  developed.  In  this  model, 
mass  transport  within  the  different  porous  structures  of  the 
DMFC  is  described  using  the  generalised  Maxwell-Stefan 
equations.  For  the  membrane  an  activity  model  based  on 
the  Flory-Huggins  approach  is  used  accounting  for  swelling 
phenomena,  related  non-idealities  and  phase  equilibria  at 
the  boundary  between  membrane  material  and  catalyst 
layer  pores.  The  model  yields  good  approximations  to 
experimental  data  with  respect  to  mass  transport  (crossover) 
and  also  reasonable  results  with  respect  to  steady-state 
current  voltage  characteristics.  It  has  to  be  pointed  out, 
that  all  simulations  were  carried  out  with  one  single  set  of 
parameters. 

The  most  significant  deviations  between  simulated  and  ex¬ 
perimental  crossover  fluxes  occur  for  high  current  densities, 
in  the  limiting  current  regime.  To  get  more  realistic  simulation 
results  in  this  respect,  two-phase  flow  in  the  anode  (carbon 
dioxide  bubble  formation)  and  cathode  (condensation  of  wa¬ 
ter  =  cathode  flooding)  pore  structures  has  to  be  accounted 
for  in  an  improved  model. 

Also  significant  deviations  exist  for  the  current-voltage 
characteristics.  This  can  be  attributed  to  the  use  of  simple 
Butler- Volmer  rate  equations  for  both  electrode  reactions. 
As  the  real  reaction  mechanisms  are  known  to  be  com¬ 
plex  reaction  networks  with  several  intermediates,  side 
and  parallel  reactions  and  coupled  adsorption/desorption 
phenomena,  a  Butler- Volmer  approach  means  a  significant 
simplification,  and  thus  can  lead  to  less  realistic  simulation 
results,  especially  in  the  kinetically  dominated  region  of 
the  current-voltage  curves  (i.e.  for  low  current  densities) 
(see  Tables  4-1 1). 

Summing  up,  further  model  refinement  and  analysis  will 
address  the  following  issues: 

•  Coexistence  of  two  phases  (gas  and  liquid)  in  the  anode 

and  the  cathode  diffusion  and  catalyst  layers. 


•  Influence  of  pressure  differences  between  anode  and 
cathode  on  the  water  and  methanol  transport  through  the 
membrane:  can  the  respective  term  for  pressure-driven 
convection  in  the  Maxwell-Stefan  equations  really  be 
neglected?  Several  experimental  studies  (also  own,  yet 
unpublished,  results)  have  shown  an  increase  in  cell 
performance  when  on  the  cathode  side  of  a  DMFC  a 
higher  pressure  is  applied  than  on  the  anode  side.  Is  this 
due  to  decreased  methanol  crossover,  or  due  to  the  higher 
oxygen  partial  pressure  at  the  cathode? 

•  Is  osmotic  pressure  (which  has  not  been  discussed  in 
this  paper  at  all)  intrinsically  accounted  for  using  the 
Flory-Huggins  activity  model? 

•  More  realistic  models  for  the  anode  and  cathode  re¬ 
action  kinetics  accounting  for  reaction  intermediates, 
adsorption/desorption  phenomena  etc. 


Appendix  A.  Physical  properties  of  pure  substances 

A.l.  Densities 


The  densities  of  all  liquid  and  solid  materials  are  assumed 
to  be  independent  of  temperature  and  pressure.  They  are  col¬ 
lected  in  Table  4.  All  gases  are  assumed  to  be  ideal,  therefore 
the  density  of  dry  air  can  be  calculated  using  the  ideal  gas 
law: 


Pair(/7,  T )  = 


M-dii 


p 


R-auT 


(A.l) 


with  the  specific  gas  constant  of  air  R a;r  =  287.221kg-1 
[31],  temperature  T  in  (K)  and  pressure  p  in  (Pa). 

The  density  of  the  mixed  anode  catalyst  can  be  calculated 
from  the  mass  fractions  w  of  both  metals  and  their  densities: 


PPtRu  =  utptppt  +  wruPru  =  0.66  ■  21400  +  0.34  •  12400 
=  18300  kg  m-3  (A. 2) 


A.  2.  Heat  capacities 

Literature  data  for  pure  substances  are  given  in  Table  5. 
For  air  as  a  standard  mixture  of  mainly  nitrogen  and  oxygen. 


Table  4 
Mass  densities 


Component  j 

Density  pj  (kgm  3) 

Liquid  water,  H2O  (1) 

997  [27d] 

Liquid  methanol,  CH3OH  (1) 

791  [27d] 

Carbon/graphite  (base  material 
of  TORAY™  carbon  paper) 

2000 

Teflon™,  PTFE 

2190  [28] 

Dry  Nafion™ 

1970  [30] 

Platinum,  Pt 

21400  [29] 

Ruthenium,  Ru 

12400  [29] 
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Table  5 

Heat  capacities 


Component  j 

Heat  capacity 

CpJ  (Jmol-1  K-1) 

CpJ  (Jkg-’K-1) 

Liquid  water,  H2O  (1)  [27b] 

75.29  [27b] 

4183  [27b] 

Carbon/graphite  (material  of 
TORAY™  carbon  paper) 

8.23  [32] 

685  [32] 

Teflon™,  PTFE 

- 

1010  [28] 

Platinum,  Pt 

- 

130  [29] 

Ruthenium,  Ru 

- 

238  [29] 

Oxygen,  02(g) 

29.36  [27b] 

- 

Nitrogen,  N2(g) 

29.13  [27b] 

- 

Water  vapour,  H20(g) 

33.58  [27b] 

- 

Carbon  dioxide  gas,  CC>2(g) 

37.11  [27b] 

- 

a  mean  heat  capacity  can  be  assumed: 

Cp. ^  =  0.79  •  Cp,n2  +0.21  •  Cp,o2  =  29.18  Jmol-1  KT1. 

(A. 3) 

The  mass-based  value  can  be  obtained  by  accounting  for 
the  molar  masses  of  oxygen  and  nitrogen: 

Cp  air  =  0.79  ■  £^+0.21  ■  =  10151kg-1  KT1. 

P '  MN2  Mo2 

(A. 4) 

A.  3.  Thermal  conductivities 

For  liquid  water  and  air  data  for  different  temperatures 
are  given  in  the  literature  (e.g.  for  water:  [33a],  for  air: 
[33b]).  The  thermal  conductivity  is  showing  a  nearly  linear 
increase  with  temperature.  Linear  regressions  yield  the  fol¬ 
lowing  simple  expressions  (with  temperature  T in  (K)  and  kj 
in  (W  m_1  K-1)): 

AH2o(i)  =  0.341  +  9.26  x  10“4  •  T  (A.5) 

A. air  =  0.0034  +  7.6  x  10“5  •  T  (A.6) 

For  wet  Nafion™,  in  the  literature  a  value  of 
A-wet Nation  =  0.43  W  m^K-1  (A.7) 

is  reported  [23]. 

A. 4.  Specific  enthalpies 

The  specific  enthalpies  hj  (Jmol-'K-1)  of  the  fluid 
components  are  calculated  from  the  standard  enthalpies 


Table  6 

Standard  enthalpies  of  formation  [27b] 


Component  j 

Ar  H?  (kJmor1) 

Liquid  water,  H2O  (1) 

-285.83 

Oxygen,  02(g) 

0 

Nitrogen,  N2(g) 

0 

Water  vapour,  H2O  (g) 

-241.82 

Carbon  dioxide  gas,  CO2  (g) 

-393.51 

of  formation,  Ap//J  (Jmol  1  K  1 )  at  standard  tempera¬ 
ture  Te  —  298.15  K,  and  the  specific  heat  capacities  Cpj 
(J  mol-1  K_1): 

hj  =  Af h]  +  CpJ(T  -  Te).  (A. 8) 

The  values  for  the  heat  capacities  are  given  in  Section 
A. 2.  The  standard  enthalpies  of  formation  are  presented 
in  Table  6. 

For  liquid  methanol,  no  standard  enthalpy  of  formation 
was  found,  but  absolute  values  of  the  specific  enthalpy  for 
different  temperatures  [33c].  A  linear  regression  and  conver¬ 
sion  from  mass  to  molar  basis  yields  the  expression 

^ch3oh(7)  =  -3726  +  48.  ST  (A.9) 

with  /icH3OH(0  in  (J  mol- 1 )  and  T  in  (K). 

A.5.  Viscosities 

According  to  [34a]  the  viscosity  of  pure  liquids  in  (Pa  s) 
can  be  calculated  from  expressions  of  the  type 

r,f  =  KT3  exp  ^A;  +  ~  +  CjT  +  DjT2)  (A.10) 

with  temperature  T  in  (K).  Table  7  presents  the  values  of  the 
parameters  for  water,  methanol  and  carbon  dioxide  as  well 
as  the  temperature  range,  for  which  they  are  valid. 

At  pressures  well  below  10  bars  (106  Pa),  the  viscosity  of 
most  gases  and  gas  mixtures  is  nearly  independent  of  the 
pressure,  but  only  a  function  of  temperature.  For  air  values 
for  different  temperatures  (at  a  pressure  of  1  bars=  105  Pa) 
are  available  in  the  literature.  These  data  show  a  nearly  lin¬ 
ear  dependence  between  viscosity  and  temperature.  A  linear 
regression  yields  the  expression 

=  (4.65  +  0.0464  ■  T)  ■  10“6  (A.ll) 

where  the  temperature  T  is  in  (K)  and  the  viscosity  results  in 
(Pa  s).  The  relative  error  compared  to  the  literature  values  is 
below  0.4%. 


Table  7 

Parameters  for  calculation  of  liquid  viscosities  [34a] 


Component  j 

T  (°C) 

Bj 

Cj 

Dj 

Water,  H20(1) 

0  to  +370 

-24.700 

4209 

0.04527 

-3.376  x  10 

Methanol,  CH3OH(l) 

-40  to  +239 

-39.350 

4826 

0.10910 

-1.127  x  10 

Carbon  dioxide,  CO2O) 

—56  to  +30 

-3.097 

48.86 

0.02381 

-7.840  x  10 
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Table  8 

Antoine  equation  parameters  for  calculation  of  vapour  pressures  [32] 


Temperature  range 

Bj 

Cj 

Water,  H20 

304-333  K 

5.20389 

1733.926 

-39.485 

334-363  K 

5.07680 

1659.793 

-45.854 

Methanol,  CH3OH 

288-356.83  K 

5.20409 

1581.341 

-33.500 

353.4-512.63  K 

5.15853 

1569.613 

-34.846 

A.6.  Vapour  pressures 


The  speed  of  sound  in  water  in  the  respective  temper¬ 
ature  and  pressure  range  is  nearly  constant  at  a  value  of 
f  sound, H20  =  1500  ms-1. 

A.  9.  Diffusion  coefficients 

Two  types  of  binary  diffusion  coefficients  are  necessary 
for  the  modeling  of  mass  transfer  using  the  Maxwell-Stefan 
approach:  One  for  each  pair  of  the  mobile  species  and  one  for 
each  mobile  species’  interaction  with  the  wall  of  the  porous 
structure.  In  general,  all  binary  diffusion  coefficients  depend 
on  the  temperature  and  the  overall  composition. 


The  vapour  pressures  of  water  and  methanol  can  be  cal¬ 
culated  using  the  Antoine  equation  [32] 

log10  (psatj  (bar))  =  Aj  -  (r(K^+c,  (A-12) 

The  parameters  A B t ,  C j  are  given  in  Table  8. 

A.7.  Liquid  molar  volumes 

For  water  and  methanol,  literature  values  are  available 
[34b]: 

Vh2o  =  18.7  x  10-6m3mol-1  (A.13) 

Vch3oh  =  42.5  x  10-6  m3  mol- 1 .  (A.  14) 

The  value  for  carbon  dioxide  can  be  calculated  using  the 
method  of  Schroeder  [34c]  as 


A.9.1.  Diffusion  coefficients  in  the  polymer  electrolyte 
membrane  (M) 

For  Nafion™,  extensive  studies  have  been  made  to  de¬ 
termine  those  diffusion  coefficients  from  experimental  data. 
In  [23]  such  expressions  are  presented.  But  as  most  re¬ 
search  in  polymer  electrolyte  membrane  fuel  cells  is  focused 
on  hydrogen-consuming  cells,  these  expressions  are  not  ac¬ 
counting  for  methanol  but  just  water  within  the  ionomer 
pores.  The  expressions  are  of  an  Arrhenius  type  to  account 
for  the  temperature  influence.  As  only  composition  influence, 
the  relative  water  content  A M  as  ratio  between  water  and 
sulfonic  acid  groups’  mole  fraction  is  accounted  for.  As  the 
sulfonic  acid  groups  at  the  pore  walls  are  not  balanced,  it  is 
more  convenient  to  use  the  proton  mole  fraction  instead  (elec¬ 
troneutrality).  The  general  expression  for  the  binary  diffusion 
coefficients  in  (m2  s-1)  according  to  [23]  is: 


Vco2  =  (3  +  2)  •  7  cm3  mol  1  =  35  x  10  6  m3mol  1 . 

(A.  15) 


Df/(L  AM)  =  D 6UAM  exp 


pA 

fl 

j_y 

R 

Te  ) 

(A.  19) 


Applying  the  same  method,  a  proton  has  a  molar  volume 
of 

Vjj+  =  7  x  10-6m3  mol-1.  (A.16) 

A.  8.  Speed  of  sound 

For  ideal  gases  the  speed  of  sound  is  related  to  the  ratio 
of  the  specific  heat  capacities  at  =  Cp/Cv,  the  temperature  T 
in  (K),  the  molar  mass  Mj  in  (kg/mol)  and  the  universal  gas 
constant  R  =  8.314  J/(mol  K)  [33e]: 


KjKl 

KsoundJ  =  W  ^  ■  (A.  17) 

For  air,  the  heat  capacity  ratio  in  the  important  temperature 
range  (300, . . . ,  400  K)  and  pressure  range  (1.5  x  105  Pa)  is 
nearly  constant  at  a  value  of  1.4  [33f]. 

The  molar  mass  of  air  can  be  approximated  assuming  air 
consisting  only  of  oxygen  and  nitrogen  with  the  respective 
mole  fractions  yo,  =  0.21  and  Tn->  =  0.79: 

Muir  ~  yo2  +  vn2  Mn2  =  28 .84  g  mol- 1 .  (A.  1 8) 


with  the  reference  temperature  Te  =  298  K.  The  values  of  the 
standard  diffusion  coefficients  and  the  activation  energies  are 
given  in  Table  2. 

As  in  the  case  of  the  DMFC  also  methanol  is  present  within 
the  membrane  pores,  three  more  binary  diffusion  coefficients 
are  necessary  for  the  pairs  (CH3OH/H2O),  (CH30H/H+)  and 
(CH3OH/M).  As  for  the  latter  two  no  literature  data  are  avail¬ 
able,  they  have  to  be  estimated  and  fitted.  As  the  methanol 
concentration  is  very  low  compared  to  that  of  water,  an  er¬ 
ror  in  these  two  diffusion  coefficients  should  have  only  a 
weak  influence  on  the  diffusion  of  water  and  protons.  But 
the  diffusion  of  methanol  is  of  course  severely  depending 
on  these  three  values.  Nonetheless,  as  methanol  is  in  many 
respects  not  so  different  from  water  (highly  polar,  small  com¬ 
pact  molecule)  as  a  first  approach  it  is  assumed  that  methanol 
has  the  same  diffusion  properties  as  water,  consequently  the 
same  parameters  are  used  as  starting  point  for  the  fitting  pro¬ 
cedure.  They  are  shown  in  Table  2. 

The  remaining  binary  diffusion  coefficient  for  the  pair 
(H2O/CH3OH)  is  calculated  assuming  a  free  solution  of 
methanol  in  water  at  infinite  dilution  (correlation  of  Hayduk 
and  Minhas  for  solutes  in  aqueous  solutions  [34d],  diffusion 
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coefficient  in  (m2  s  1)): 
tV'HiOH.lhO  ~  £>Oi3OH,H20 

=  1.25  x  10-12(ycH3oH  -  0.292)7’1-52(??™or  (A.20) 

with  the  exponent 
_  9.58 

e*  =  = - 1.12,  (A. 21) 

Vch3oh 

the  molar  volume  of  methanol  Vch3oh  =  42.5  cm3  mol-1 
and  the  viscosity  of  pure  water  rj^O  111  (c?  =  103  Pa  s). 

A.9.2.  Diffusion  coefficients  in  the  anode  diffusion  layer 
(AD) 

The  binary  diffusion  coefficients  of  the  mobile  species 
in  the  anode  diffusion  layer  (AD)  are  calculated  using  the 
Tyn-Calus  method  [34e]for  diffusion  coefficients  in  liquid 
solutions  at  infinite  dilution  (in  (m2  s-1)): 

=  *-93  *  >o-'2  (|l)  (A-22) 

Here  component  i  is  the  solute  and  j  is  the  solvent.  The 
molar  volumes  Vj  are  in  (cm3  mol-1),  the  viscosities  tA,s  in 
(cP  =  103  Pas)  and  the  temperature  T  in  (K).  P,  and  Pj  are 
so-called  parachors,  which  are  related  to  the  liquid  surface 
tension,  but  can  also  be  estimated  from  a  groups  contribution 
method  developed  by  Quayle  [34e].  For  water,  methanol  and 
carbon  dioxide,  this  method  leads  to  parachor  values  of 

Ph2o  =  2  •  PH  +  P-O-  =  2  x  15.5  +  20 

=  51cm3ga25s-0'5,  (A. 23) 


Pch3oh  =  Pch3  +  P-oh  =  55.5  +  29.8 

=  85.3cm3g°'25s-0-5,  (A. 24) 

Pco2  =  Pc  +  2  •  Po  =  9  +  2  X  20  =  49 cm3  g025  s-05. 

(A. 25) 

According  to  the  literature,  if  water  is  the  solute,  the  para¬ 
chor  and  molar  volume  values  of  water  shall  be  doubled  (wa¬ 
ter  is  treated  as  a  dimer). 

With  the  help  of  Eq.  (A.22),  all  six  diffusion  coefficients 
at  infinite  dilution  for  the  three  species  (water,  methanol,  car¬ 
bon  dioxide)  can  be  calculated.  To  get  the  necessary  three 
binary  diffusion  coefficients,  each  pair  of  the  former  six  val¬ 
ues  belonging  to  the  same  two  species  are  combined  using 
the  method  of  Vignes  [34f] : 

Dy  «  Dr,  =  « Vi gnes  [ ( Df°  )XJ  ( Dfjf  )■*'].  (A. 26) 

Both  mole  fractions  x,-  and  xj  are  set  to  0.5,  the  thermo¬ 
dynamic  factor  avignes  is  assumed  to  be  1  (ideal  mixing  be¬ 
haviour  of  the  two  species),  leaving  the  expression 

D u  =  (D’-j  DjP)0'5.  (A. 27) 


As  the  mass  transport  takes  place  within  a  porous  ma¬ 
trix,  effective  diffusion  coefficients  are  needed.  To  convert 
the  gained  values  into  effective  coefficients,  it  has  to  be  ac¬ 
counted  for  the  morphology  of  the  solid  matrix  represented 
by  the  porosity  e  and  the  tortuosity  coefficient  r  [18].  (Note 
that  there  is  a  mistake  in  the  mentioned  reference:  there  both 
diffusion  coefficient’s  indices  have  been  confused.) 

DfF  =  -Dy.  (A. 28) 

J  r 

To  describe  the  ratio  between  tortuosity  coefficient  and 
porosity,  many  approximations  exist.  According  to  [18],one 
of  the  most  commonly  used  is  based  on  the  approximation 
that  the  tortuosity  is  only  a  function  of  the  porosity  (and  not 
of  the  size  of  the  mobile  species): 

r  =  e-L5.  (A. 29) 


A.9.3.  Diffusion  coefficients  in  the  cathode  diffusion 
layer  ( CD) 

In  the  cathode  diffusion  layer  it  is  assumed  that  all  mobile 
species  (oxygen,  nitrogen,  water  vapour  and  carbon  dioxide) 
are  ideal  gases.  The  diffusion  coefficients  describing  the  in¬ 
fluence  of  the  pore  wall  are  calculated  using  the  Knudsen 
equation  (according  to  [18]), 


eff  _  e  ^pore  I  $RT 

'/M  ~  T~\JnWi’ 


(A. 30) 


although  here,  strictly  speaking,  no  Knudsen  diffusion  takes 
place  due  to  the  big  mean  pore  diameter  in  the  carbon  pa¬ 
per  which  is  some  orders  of  magnitude  bigger  than  the  gas 
molecule  diameters.  It  turns  out  that  the  coefficients  calcu¬ 
lated  from  Eq.  (A. 30)  are  of  such  an  order  of  magnitude,  that 
in  the  end  there  is  no  significant  influence  of  the  wall  friction 
on  diffusion. 

The  pair  diffusion  coefficients  of  the  mobile  species  are 
calculated  according  to  the  method  of  Fuller,  Schettler  and 
Giddings  in  free  gas  phase  [33d](in  (m2  s-1)): 

IQ-1  T115  J  (Mi  +  Mj)/ Mi  Mj 

DU  *  DU  = - ^3 - ;  1/3  — ■  (A. 31) 

p[(£  V*)J/3  +  (£  V*)/3]2 

Here  the  temperature  7’ is  in  (K),  M  are  the  molecular  weights 
in  (g/mol),  p  is  the  pressure  in  (atm=  105  Pa)  and  the  sum 
terms  are  sums  of  atomic  diffusion  volumes,  which  are  tab¬ 
ulated  for  many  simple  gases  (Table  9).  To  get  effective  dif¬ 
fusion  coefficients,  also  here  Eq.  (A. 28)  is  used. 

In  the  simulations,  the  presented  diffusion  coefficients 
were  used  as  starting  values.  In  order  to  get  a  better  fit  to 
experimental  data,  some  of  these  values  have  been  slightly 
adjusted,  as  presented  in  Table  2. 

Table  9 


Parameter  values  for  Eq.  (A. 3 1)  (atomic  diffusion  volumes  taken  from  [33d]) 


o2 

n2 

h2o 

C02 

Mj  (g  mol  1 ) 

32 

28 

18 

44 

CT,  V*)j  (cm3  mol-1) 

16.6 

17.9 

12.7 

26.9 
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Table  10 

Physical  properties  of  TGP-H-060  carbon  paper  [35] 


Thickness  (|xm) 

Electrical  resistivity  (m£2  cm) 

Thermal  conductivity  (W  m  1  K  1 ) 

Porosity  (-) 

Bulk  density 
(g  cm”3) 

Through  plane  In  plane 

Through  plane  (20  °C)  In  plane  (20  °C) 

in  plane  (100  °C) 

190  (170)a 

80  5.8 

1.7  21 

23  0.78 

0.44 

aOwn  measurement  used  for  further  calculations. 


Appendix  B.  Porosities  and  volume  fractions 

The  major  physical  properties  of  TORAY™  paper  TGP- 
H-060  (according  to  data  sheet  supplied  by  Toray  Deutsch¬ 
land  GmbH)  are  given  in  Table  10.  As  the  carbon  paper  is 
PTFE  treated  prior  to  use  in  the  DMFC,  the  real  porosity  has 
to  be  calculated  as  a  function  of  the  PTFE  content.  The  to¬ 
tal  volume  of  the  carbon  paper  is  the  sum  of  the  carbon,  the 
PTFE  and  the  free  pore  volume: 


and  PTFE,  £ptfe  =  2.19g/cm3  and  gcarbon  =  2.0  g/cm3  the 
porosity  of  the  PTFE-treated  TORAY™  paper  results  as 

eporef_treatedToiay  =  0.71.  (B.9) 

The  PTFE  volume  fraction  is  then  finally 


PTFE-treated  Toray  _  ,  PTFE-treated  Toray 

fcPTFE  —  1  b  pores 

PTFE-treated  Toray 
-  ^carbon  =  °-07 


fB.  10) 


Vtot  =  Fcarbon  +  VpTFE  +  Ppores 


(B.l) 


The  carbon  volume  can  be  expressed  by  the  porosity  of 
the  untreated  carbon  paper  (see  Table  10): 


^untreated  Toray  \  t  / 
epores  >  Hot 


fcarbon  —  (1 

The  carbon  volume  fraction  is  consequently 


(B.2) 


untreated  Toray  PTFE— treated  Toray 

^carbon  ^carbon 

_  1  _  untreated  Toray  _ n  99 

1  Spores  u.zz. 

The  PTFE  mass  fraction 
m  PTFE 


WPTFE  = 

^carbon  +  /«PTFE 

can  be  rearranged  to  get  the  PTFE  mass 
WPTFE 

"I  PTFE  =  - - m  carbon 

I  —  Ut  ptfe 


(B.3) 


(B.4) 


(B.5) 


Substituting  masses  by  bulk  densities  times  volumes  and 
rearranging  we  get  the  PTFE  volume 


WPTFE  Pcarbon 

PPTFE  —  Fcarbon- 


(B.6) 


1  —  IP  PTFE  PPTFE 
Substituting  the  carbon  volume  using  Eq.  (B.2)  one  gets 

WPTFE  Pcarbon  , ,  .untreated  Toray sy 


PPTFE  = 


_  untreai 

Spores 

1  —  W  PTFE  PPTFE 


(B.7) 


The  resulting  porosity  (i.e.  the  pore  volume  fraction)  of 
the  PTFE-treated  material  is  then 

PTFE-treated  Toray 
spores 


_  spores  _  untreated  Toray 

T  /  spores 

Uot 

_  WPTFE  Pcarbon,.  _  untreated  Toray-. 
i  ^  spores  1‘ 

1  —  W PTFE  PPTFE 


(B.8) 


With  the  parameters  given  in  Table  11,  the  typical  PTFE 
mass  content  u>ptfe  =  0.25  and  the  bulk  densities  of  carbon 


Also  the  porosity  of  the  catalyst  layers  has  to  be  deter¬ 
mined.  It  can  be  calculated  as  the  sum  of  the  volumes  of 
all  dry  materials  within  the  catalyst  layer  (i.e.  catalyst  and 
NAFION™)  divided  by  the  total  volume  of  the  catalyst  layer. 
The  catalyst  and  NAFION™  volumes  are  calculated  from  the 
masses  and  the  densities: 


total  volume  —  catalyst  volume  —  NAFION  volume 
total  volume 


(B.l  1) 


One  ends  up  with  an  equation  using  the  catalyst  and 
NAFION™  loadings  wcat  and  mj  nafion,  and  the  catalyst  layer 
thickness  c/caLlayer: 


l/AS(Vt“t'layer  -  (/Hcat/ Peat)  -  (/HNAFION/PNAFION)) 


/^cat-layer  _  (wcat/pcat)  -  (  WNAFIOn/PNAFIOn) 

^cat.layer 


(B.12) 


Table  1 1 

Physical  properties  of  catalyst  layers  and  calculation  of  porosities 


Property 

AC 

cc 

Thickness  (rfcat.layer) 

35  (Jtm 

35  (Jim 

Catalyst  loading  (wc at) 

5  mg  cm-2 

5  mg  cm-2 

Catalyst  density  (yOcat) 

0.66  x  21.4  +  0.34  x 

21.4 gem”3  (Pt) 

NAFION™  loading 

12.4  =  18.3  gem”3 
(Pt:Ru  mass  ratio 

66:34) 

0.15  lOcat  = 

0.10  Wcat  = 

(u/nafion) 

0.75mg  cm-2 

0.5  mg  cm-2 

NAFION™  density 

1.97  gem”3 

1.97  gem”3 

(Pnafion) 

Calculation  of  volume  fractions 

Catalyst 

0.08 

0.07 

(solid  metal) 

NAFION™ 

0.11 

0.07 

Free  pore 

0.81 

0.86 

space  =  porosity  (e) 
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For  the  applied  MEAs  the  resulting  calculations  are  shown 
in  Table  1 1 . 


Appendix  C.  Effective  thermal  conductivities 

From  Toray  Corporation  a  value  for  the  through-plane 
thermal  conductivity  is  available  (see  Table  10).  But  this  value 
is  only  valid  for  the  pores  being  filled  with  air,  at  20  °C.  There¬ 
fore,  this  value  is  not  applicable  for  the  anode  diffusion  layer 
(AD)  where  the  pores  are  filled  with  water  (and  traces  of 
methanol  and  carbon  dioxide).  Also,  as  the  thermal  conduc¬ 
tivity  of  air  is  temperature  dependent,  using  this  value  for  the 
cathode  diffusion  layer  seems  problematic.  Therefore,  first 
the  thermal  conductivity  of  the  carbon  fibres  alone  shall  be 
calculated.  Generally,  assuming  parallel  thermal  conduction 
through  all  present  materials,  the  effective  thermal  conduc¬ 
tivity  can  be  calculated  from 

*eff  =  E£;V  (c-d 

j 

Using  this  expression  for  TORAY™  paper  at  20  °C  we  get 
the  (assumed  temperature-independent)  thermal  conductivity 
of  TORAY™  paper  in  the  vacuum  (i.e.  the  contribution  from 
the  carbon  material  only): 

,  eff ,  20°  C  untreated  Toray ,  20°  C 

,  _  '''Toray  +  air  Spores  Ciir 

Toray  untreated  Toray 

‘  Spores 

=  7.63  W  m-1  K-1  (C.2) 


Appendix  D.  Volumetric  effective  heat  capacities 

In  the  energy  balances  of  the  simulation  model,  the  vol¬ 
umetric  overall  heat  capacities  in  all  control  volumes  are 
needed.  These  are  calculated  from  the  densities  pj  (kg  m-3), 
the  mass-based  heat  capacities  C p.  j  (J  kg- 1  K_1)  and  volume 
fractions  sj  (-)  of  all  present  materials  and  reactants: 

(pCp)  —  'y  'sjPjCpj.  (D.l) 

j 

The  pores  in  the  anode  diffusion  layer  (AD)  are  assumed  to 
be  filled  with  liquid  water  and  small  amounts  of  methanol  and 
carbon  dioxide.  As  the  mole  fractions  of  the  latter  two  are  each 
well  below  0.05,  their  influence  on  the  overall  heat  capacity 
is  neglected.  It  is  only  accounted  for  the  carbon  material  and 
the  PTFE  forming  the  solid  porous  structure  as  well  as  the 
water  inside  the  pores: 

,"pr\AD  _  PTFE— treated  Toray  ^ 

xP'-'p)  —  ^carbon  Pcarbon^/?,  carbon 

PTFE— treated  Toray  ^ 

+  £pTFE  PPTFE  C  PTFE 

+  £pores>trCated  T°ray  PH20(1)  C  ;p  h20(1)  •  (D.2) 

The  identical  equation  is  applied  to  the  cathode  diffusion 
layer  (CD).  The  only  difference  is  that  the  pores  are  filled 
with  air  instead  of  liquid  water. 

As  in  the  anode  diffusion  layer  (AD),  all  pores  in  the  anode 
catalyst  layer  (AC)  are  assumed  to  be  filled  with  pure  water. 
Methanol  and  carbon  dioxide  are  again  neglected.  Here  the 
solid  matrix  is  formed  from  the  catalyst  particles  and  ionomer 
(Nafion™): 


The  thermal  conductivity  of  air  at  20  °C  is  calculated  from 
Eq.  (A. 6).  The  other  values  are  taken  from  Table  11. 

Using  the  value  calculated  in  Eq.  (C.2)  and  the  volume 
fractions  of  carbon  fibres  and  PTFE  as  well  as  the  porosity 
(see  Appendix  B),  the  effective  thermal  conductivity  can  be 
calculated  accounting  for  the  material  filling  the  pores  and 
for  the  temperature. 

In  the  anode  diffusion  layer  (AD)  the  pores  are  assumed 
to  be  filled  with  a  liquid  mixture,  which  mainly  consists  of 
water.  For  simplicity,  the  influence  of  methanol  and  carbon 
dioxide  is  neglected.  The  resulting  expression  is 


^eff.AD 


PTFE— treated  Toray .  PTFE— treated  Toray . 

=  ^carbon  AToray  +  £PTFE  APTFE 


+  £poTref“treatedTOray^H2o(r)  (C.3) 


For  the  thermal  conductivity  of  water,  the  literature  value 
Eq.  (A. 5)  is  used. 

In  the  cathode  diffusion  layer  (CD)  the  pores  are  assumed 
to  be  filled  with  air.  Eq.  (A.6)  is  used  for  the  thermal  conduc¬ 
tivity  of  air. 


^eff.CD 


PTFE— treated  Toray PTFE— treated  Toray , 

=  ^carbon  AToray  +  ^PTFE  APTFE 


PTFE— treated  Tor 
'  pores 


^ahCO 


(C.4) 


(pCp)AC  —  £ptRu(pCp)ptRu  +  £Nafion(P^'/>)Nafion 

+  epOTes/°H20(/)Cp,H20(/)  (D.3) 

with  the  heat  capacity  of  the  PtRu  catalyst  calculated  from 
the  mass  fractions  of  both  metals: 

(pCp)ptRu  =  0.66pptCp,pt  +  0-34pRuCp.Ru 

=  2.84  x  106Jm“3K_1  (D.4) 

In  the  literature  [23]  for  wet  NAFION™  one  finds  a  vol¬ 
umetric  heat  capacity  of 

(pCp) Nation  =  2.4  x  107J  m^K”1 .  (D.5) 

For  the  cathode  catalyst  layer  (CC)  the  same  reason¬ 
ing  is  valid  as  for  the  anode  catalyst  layer  (AC),  except 
that  the  pores  are  filled  with  air,  and  the  catalyst  is  pure 
platinum. 


Appendix  E.  Concentration  measures  within  the 
PEM 

To  convert  the  concentration  measures  used  in  the  PEM 
model,  the  concentration  with  respect  to  the  total  volume 
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including  the  polymer  backbone  c,  (mol  m-3)  is  useful.  It 
can  be  calculated  from  the  molar  density  Nf  in  the  spatially 
discretised  model  as 


'•M  = 

J,k 


A41 


(E.l) 


with  j  =  H2O,  CH3OH,  H+,  where  the  thickness  (m) 
of  the  respective  control  volume  k  is  calculated  from 

r/M'dry 


AZiM  = 


ncvJ 


M 


£ 


(E.2) 


This  thickness  is  also  needed  for  the  various  flux  calcu¬ 
lations  to  calculate  the  distances  between  the  centres  of  the 
control  volumes.  From  Cj,  all  other  concentration  measures 
can  be  easily  calculated: 


volume  fractions  :  e^1  =  (f  V j . 


-M 


concentrations  within  the  pores  :  — 


CM  ’ 
"pores 


with 

eM 

pores 


E  M  M  ,  M  ,  „M 

Ej  —  £H+  +  eH20  ■+■  fcCH3OH' 


(E.3) 

(E.4) 

(E.5) 
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